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Abstract: We present the two first terms in the threshold expansion of Higgs production 
partonic cross-sections at hadron colhders for processes with three partons in the final 
state. These are contributions to the inclusive Higgs cross-section in gluon fusion at N3LO. 
We have developed a new technique for the expansion of the squared matrix-elements 
around the soft limit and for the reduction of the required phase-space integrals to only 
ten single-scale master integrals. We compute the master integrals building upon modern 
techniques for the integration of multidimensional integrals in dimensional regularization. 
Our results constitute an important step towards a systematic computation of the Higgs 
boson cross-section as an expansion around the threshold limit. 



Keywords: QCD, NLO, NNLO, N3LO, LHC, Tevatron. 



1. Introduction 



The first years of experiments at the Large Hadron Cohider (LHC) have probed decisively 
the physics of electroweak symmetry breaking, demonstrating a great understanding of the 
theory at this energy scale. The success of particle theory must be attributed not only to 
the Standard Model (SM), but also to remarkable progress in perturbative calculations. 
Results in this field of research are outstanding and have resulted in very accurate com- 
parisons with data. Indeed, partonic cross-sections for a large variety of processes are now 
routinely computed at next-to-leading order (NLO) in the strong coupling constant ag ^. 
In addition, modern Monte-Carlo programs have reached a high level of sophistication with 
the development of methods to match and merge parton showers with fixed-order calcula- 
tions^. Finally, basic processes at the LHC with a small final-state multiplicity can now be 
computed fully differentially at next-to-next-to-leading order (NNLO)^. 

A classic example where perturbative calculations have been crucial for the inter- 
pretation of LHC data is the production of a Higgs boson. Perturbative corrections for 
Higgs boson signal processes are generically large, especially so for production via gluon 
fusion [49-55]. For this process, the total and fully differential cross-sections are known 
through NNLO with a typical uncertainty of about 10% due to variations of the factor- 
ization and renormalization scales [46,56]. Indeed, without the knowledge of the NNLO 
corrections, predictions at LO or NLO would be assigned theory uncertainties which are 
larger than the already achieved experimental uncertainties. With future LHC data, a 
comparison with theory at a level of precision of a few percent will be revealing the fine 
details of the mechanism of electroweak symmetry breaking and potentially uncovering 
new fundamental laws of physics. 

Besides being of importance to Higgs physics, NNLO predictions have promoted the 
Drell-Yan process to a standard candle for studies at hadron colliders, where observables 
in this process arc computed with a typical precision of better than 2% as indicated by 
scale variations [57, 58] . Measurements of the mass of the W boson and the weak mixing 
angle from Drell-Yan data are one of the most stringent constraints on the Standard Model. 
Furthermore, Drell-Yan production data is an essential input for the extraction of parton 
distribution functions from hadron collider processes. Last but not least, the clean detector 
signatures of Drell-Yan events, combined with the excellent theoretical predictions of their 
rates, allow for a precise determination of the luminosity. Indeed, luminosity determination 
from Drell-Yan is relatively insensitive to high pile-up conditions, a fact which will be even 
more important during the high-luminosity phase of the LHC. 

We believe that it is possible to develop efficient methods for computing Higgs boson 
and Drell-Yan hadroproduction cross-sections at the next perturbative order, "N3LO" . The 
benefits of such a program are potentially very important. 

In Higgs production via gluon fusion the uncertainty due to scale variations at N3LO 
is anticipated to be half of what is found at NNLO [59]. Scale variation estimates can- 

^See, for example, refs. [1-14] 

^See, for example, refs. [15-35] 
^See, for example, refs. [36-48] 
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not be derived from first principles, but entail a certain degree of subjectiveness. With 
the knowledge of more orders in the perturbative expansion, the remaining perturbative 
uncertainty can be estimated more reliably, not only from scale variations but also from 
the progression of the series. Improvements of the Higgs boson production cross-section 
uncertainty will propagate into the determination of the couplings of the Higgs boson. 

Although existing NNLO calculations for the Drell-Yan process indicate that the cross- 
section is already known very precisely, an explicit N3LO computation will educate us 
further by challenging or verifying the traditional prescriptions for uncertainty estimates. 
Furthermore, high-luminosity runs at the LHC require the triggering of Icptons at a higher 
transverse momentum, resulting in a deterioration of the NNLO scale uncertainty [45] . An 
excellent theoretical precision can then be recovered by including the N3LO perturbative 
corrections. 

So far, there has been no N3LO calculation for hadron collider processes performed in 
the literature. In this paper, we attempt a first step towards computing inclusive hadropro- 
duction cross-sections at this perturbative order. We focus on one of the most complicated 
N3LO contributions to the inclusive Higgs boson cross-section in the gluon fusion channel, 
namely, the contributions coming from partonic cross-sections for tree-level radiative pro- 
cesses with three additional partons emitted in the final state. A direct integration of the 
corresponding matrix-elements over phase space is however tedious. To simplify the prob- 
lem, we perform a threshold expansion and devise a method to compute the coefficients 
of the expansion analytically. One of the main results of our paper is the computation of 
the first two terms of the threshold expansion for triple real corrections to inclusive Higgs 
production. 

Our expansion method builds upon reverse unitarity [55,58,60-62], a technique devel- 
oped for the calculation of the inclusive Higgs cross-section and the rapidity distribution of 
electroweak gauge bosons at NNLO. The main idea is that phase-space integrals of matrix- 
elements are dual to loop integrals in their algebraic properties: recurrence identities in 
the powers of propagators and the number of dimensions as well as differential equations 
satisfied by loop integrals also apply to phase-space integrals. This is achieved by associat- 
ing on-shell conditions or other phase-space constraints in the form of delta functions with 
differences of otherwise identical propagators with opposite infinitesimal imaginary parts, 
in the spirit of Cutkosky's rules. In this paper, we observe and exploit that the duality 
rules of reverse unitarity can be expanded in kinematic parameters. 

After performing a threshold expansion of loop integrals dual to the cross-sections for 
Higgs plus three partons processes, we apply automated reduction algorithms to reduce 
the coefficients of the expansion to master integrals. We have performed the reduction 
using existing public programs and programs developed specifically for the purpose of this 
computation. 

The reduction yields ten master integrals for the first two coefficients in the threshold 
expansion. The master integrals themselves are not specific to Higgs production and will 

appear in the threshold expansion of other hadroproduction processes, such as Drell-Yan. 
For processes with a single colorless final state our set of master integrals is complete and 
our expressions are universal for the partonic cross-sections with the same initial state at 
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leading order in the threshold expansion. 

For the computation of the soft master integrals we have employed a combination of 
methods. The master integrals are evaluated by using a parametrization of the phase space 
in terms of the energies and angles of the momenta of the soft partons. When needed, 

we perform the phase-space integrals by introducing Mellin-Barnes representations. We 
evaluate some of the McUin-Barnes master integrals in a closed form as hyper geometric 
functions or directly as a Laurent-expansion in the dimensional regulator e with infinite 
nested sums as coefficients, which naturally evaluate to multiple zeta values in all cases. 
Some master integrals cannot be evaluated easily from their Mellin-Barnes representations. 
We have developed a procedure to turn Mellin-Barnes integrals into integrals over positive 
real parameters, which are easy to expand in e whenever the integral is finite. The resulting 
parametric integrals are then evaluated by integrating out the integration variables one-by- 
one in terms of multiple polylogarithms. The last integration then again produces multiple 
zeta values in a natural way. 

We have restricted our computation of the partonic cross-sections to the first two terms 
of the threshold expansion. There are two obvious ways to extend this work in the future. 
We could continue with the computation of sublcading terms in the threshold expansion. 
This requires a more intensive but not prohibitive computer algebra and the evaluation 
of soft master integrals from topologies which contribute only to higher orders in the soft 
expansion. The soft expansion is rather fast converging for Higgs production but rather 
slow for Drell-Yan, as has been noticed at NNLO. A second possibility is to perform the 
phase-space integrations for arbitrary kinematics, reducing them to a different set of master 
integrals. The soft master integrals which we present in this article can serve as boundary 
conditions for solving the differential equations satisfied by the master integrals in arbitrary 
kinematics. 

This article is organized as follows: In Section 2 we present our general strategy based 
on reverse-unitarity to perform the threshold expansion for real-emission cross sections in 
terms of soft integrals, and we illustrate the method on some simple examples in Section 3. 
In Sections 4 and 5 we study some properties of soft integrals in general, and we show 
that there is an easy and canonical way to derive dimensional recurrence relations and 
Mellin-Barnes integral representations for generic soft integrals. The soft master integrals 
contributing to the first two terms in the threshold expansion of the leading-order partonic 
cross sections for pp ^ H + 3 partons are discussed in Section 6, and the corresponding 
results for the cross sections are presented in Section 7. Section 8 contains technical details 
about the computation of the soft master integrals, and in Section 9 we draw our conclusions 
and give an outlook for future work. The paper contains several appendices discussing phase 
parametrisations and a generic formula for the phase-space volume foi H + n partons, as 
well as a new method to derive a parametric integral representation from a Mellin-Barnes 
integral and a description of an algorithmic way to perform the analytic integration of 
certain classes of parametric integrals. 
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2. Reverse unitarity, threshold expansion and soft integrals 



We consider the production of a Higgs boson in association with j = 3. . .N massless 
partons in the final state from two massless partons i = 1, 2 in the initial state, 



1 + 2^ H + 3 + ... + N. 



(2.1) 



The inclusive cross-section for this process in dimensional regularization is given by a 
phase-space integral over the momenta qj of the final-state partons, 

a = J d<i>N-iiqH,q3,...,qN;M^;s;D)\Af i{qj},qi,q2;D). (2.2) 

We work in D = 4 — 2e dimensions and denote by qi , q2 the momenta of the initial-state 
partons and we also use the shorthand notation 

qi2 = qi+ q2, q'345 = q3 + q4: + qs, etc. 

In the following we will often drop the functional dependence on the dimension D for clarity. 
The mass of the Higgs boson is denoted by M, and we denote the (squared) center-of-mass 
energy by s = (gi -|- (?2)^. |.4p represents the squared matrix-element multiplied with the 
appropriate flux and symmetry factors. The D-dimensional phase space measure is given 

by 

d^N-iiqH,q3, ■■■,qN] M^; s- D) 



j=3 



with (5+(g^ — m?) = 5{q^ — m^)9(g°). Integrating out the momentum of the Higgs boson, 
we can rewrite eq. (2.2) as 



a = 



27r 



N 



i=3 



{[q3...N - qi2f - M^) \Af {{qj},qi,q2). (2.4) 



In this article, we restrict ourselves to the case of real-radiation matrix-elements with- 
out virtual corrections. We introduce the variables 



M2 

z = and z = 1 — z . 

s 



(2.5) 



We now rescale the momenta of all the partons. 



qi 



y/spi, ifi = l,2, 
^/s zpi , ii i = 3 . . . N , 



(2.6) 



which captures the scaling of the partonic momenta in the final-state. We emphasize that 
this is not, as yet, an approximation, but rather a convenient change of integration variables 
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which captures the correct asymptotic behavior at threshold as z ^ 0. In the following we 
assume s = 1, and we find 



a 



-{D-2){N-2)-l 



27r 



N 

n 

j=3 



(2^ 



\P3...N-Pl2f - ZpI_n) 



(2.7) 



X 1^1^ {{zpj},Pi,P2)- 



Note that the full s-dependence can easily be recovered from dimensional analysis. 

The squared matrix-element \ Af consists of a rapidly growing number of terms with N, 
yielding a correspondingly large number of phase-space integrals. The method of reverse 
unitarity, developed in refs. [55,58,60-62], allows the reduction of phase-space integrals 
to a basis of fewer master integrals by establishing a duality of phase-space and loop 
integrals, where the latter are amenable to algebraic methods [63, 64] based on integration 
by parts [65, 66]. 

Following the reverse unitarity approach, on-shell and other phase-space constraints 
are dual to propagators: 



1 1 
— ^ Disc 



2m \q'^+iO 



(2.8) 



"Cut" propagators can be differentiated in a similar way to ordinary propagators with 
respect to their momenta. 



d 



1 



2g^ 



(2.9) 



leading to identical integration-by-parts (IBP) identities for phase-space integrals as for 
their dual loop integrals. Solving the system of IBP identities for phase-space integrals 
proceeds in the same way as for loop integrals, with the exception that for cut-propagators 
we can use the simplifying identity: 



1 



-^0, Vz/ = 0,1,2,... 



(2.10) 



According to the reverse unitarity method, we find a dual forward scattering loop- 
amplitude with N — 1 cut-propagators for the real radiation contribution of eq. (2.7), 
namely. 



-(D-2)(Ar-2)-l J_ 



27r 



N 

n 

J=3 



(27r) 



D-l 



1 



\PZ...N-Pl2f - Zpl_^ 



(2.11) 



X \A\^{{zpj},puP2). 



In this article, we take one further step and expand cut-propagators and the squared 
matrix-elements around z = 1, 



A;=0 



(2.12) 
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and 



\P3...N - P12] - f^Q [pi2 • (Pl2 - 2p3...iV 

In this approximation, the cross section can be expanded into a power series in z, 



^ ^ ^(D-4)(iV-2)- 



-i^^fc^5(fe)_ (2.14) 



A;=0 

The coefficients of the power series are given by 



a'^'^ = t.(-^y f „ V {pLNf\-A\U{{Pj},PuP2), (2.15) 

J lPl2 ■ {Pl2 - ^P3...N)j c 

where d<^'fj_^ denotes the "soft" phase space measure 

n 



27r [pi2 ■ (Pi2 - 2i?3...Ar) 
— S+ipl2 - 2pu ■ P3...N) n (27r)^-i '^+^^- 




TV 

2^ 

'j' 

3=3'^" ' 



(2.16) 



The integrals which emerge after the z expansion depend trivially on one dimensionful 
parameter = If we put s = 1, the integrals are number integrals whose only functional 
dependence is through the space-time dimension D = 4 — 2e. We will refer to such integrals 
as soft (phase space) integrals, and they arc the main subject of this paper. Wc note that, 
apart from the cut Higgs boson propagator, the integrands of soft phase space integrals 
are homogeneous functions under a simultaneous rescaling of the final-state momenta. In 
addition, a soft integral can be reduced to a set of "soft" master integrals using IBP 
identities by exploiting the duality to loop integrals via reverse-unitarity. We will illustrate 
this property in the next section where we check our method on several examples. 

3. Validation of the method and examples 

In this section, we study the validity of the method described in the previous section at 
NLO and NNLO - two perturbative orders that are well studied in the literature and so 
we can compare our results readily with known results. In particular, we show that our 
method reproduces the correct results for the leading behavior of NLO and NNLO real 
emission amplitudes in the soft limit, as well as for the subleading terms in the expansion 
of the phase space volume up to N3LO and for a non-trivial double real emission master 
integral at NNLO. 

At NLO, all phase space integrals that contribute to the real emission amplitude in 
general kinematics can be reduced to the phase space volume for + 1 parton, 

e) = , \, ^/^ ~ '\ z'-^' . (3.1) 
^ ^ 2(47r)i-^ r(2-2e) ^ ^ 
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As there is only one master integral which is a monomial in z, our method trivially gives 
the correct answer at NLO. 

At NNLO all double real emission phase space integral can be reduced in general 
kinematics to a linear combination of 18 master integrals [55]. The leading contribution 
of all master integrals in the soft limit to all orders in e was computed in ref. [67], and it 
was observed that in this limit 17 master integrals are proportional to the soft limit of the 
phase space volume for H + 2 partons, 

*»<^^^' =W^'"~"W^f^'^'-'-'-^-''-''-''^ (3.2) 



where we defined 



More precisely, it was shown in ref. [67] that if Xf(2:;e) denotes the leading term in the 
soft limit of the double real emission master integrals, then we can write"^ 



Xf(f;e) =Si(f;e)$t(e), 1 < i < 17, 

Xf,(.; .) = -4 - ^''<^ - ^'>'^ - 3F.(1. 1. 1 - , 1 - 2.; 1) (.) . 

where Si(z; e) are monomials in z and rational functions of e. Using the method described 
in the previous section, we can easily explain the structure of eq. (3.4). Indeed, we observe 
that in the soft limit all the double real emission phase space integrals can be reduced 
to only two master integrals. In particular, the IBP identities in the soft limit allow us 
to express all but one of the in terms of the phase space volume, and the coefficients 
appearing in the reduction are precisely the functions Sj. In other words, in the soft limit 
all double real emission phase space integrals can be reduced to linear combinations of the 
following two soft master integrals 




d^l , (3.5) 



5 — . (3.6) 

S14S23S34 



Our method thus provides the correct leading soft behavior of the double real emission 
contribution at NNLO. We emphasize that all the diagrams in this paper represent soft 
phase space integrals, i.e., all the diagrams represent integrals with respect to the soft 



*Note that the normahzation differs sUghtly from the normalization of ref. [67], 
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phase space measure of eq. (2.16). In addition, the invariants appearing in the integrands 
of the soft integrals are defined with respect to the rescaled momenta defined in eq. (2.6), 




{npi + TjPj) , = { ,r.' \r (3.7) 



Our method does not only allow us to compute the leading soft behavior, but we can 
consistently expand around the soft limit z = 0. In the following we show that we can 
correctly reproduce the first few terms in the soft expansion of double and triple emission 
phase space volumes, as well as for the NNLO master integral Xig of refs. [55,67,68]. 

Let us start with the phase space volume for H + 2 partons in the limit where the two 
partons are soft. On the one hand, from eq. (3.2) we immediately see that $3 admits the 
expansion 



3-4e ^S^,^ (1 ~ e)n(2 - 2e)„ _^ 



n=0 (4 - 4.)„ 



S4.^'?. 1-e (l-e)(2-e)(3-2e) o .^r 
= z^-^^ a>f (e) \^l + ^-z+ ^ ^ + . 

On the other hand, using eq. (2.13) we obtain the diagrammatic expansion 



(3.8) 



^3(^;e) 



+ o(r 



(3.9) 



where the dashed lines indicate numerator factors and dots represent additional powers 
of the propagators or the numerators. The diagrams appearing in eq. (3.9) are in one-to- 
one correspondence with the terms in the expansion (3.8). Indeed, IBP reduction of the 
integrals in eq. (3.9) reveals 




:i-e)(2-6)(3-26) 
4(5 - 4e) 



3^ 



(3.10) 
(3.11) 



in perfect agreement with eq. (3.8). We checked explicitly that our method reproduces 
correctly the first ten terms of the soft expansion of the phase space volume for H + 2 
partons. 

As a second example we derive the subleading terms in the soft expansion of the 
double real emission master integral Xig. Unlike for the phase space volume, no result is 
known for Xig valid to all orders in e in general kinematics, but the integral was evaluated 
explicitly up to 0{e) in terms of harmonic polylogarithms [69] in ref. [55,67,68]. We can 
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thus compare the result of our method order by order in e to the expansion of the harmonic 
polylogarithms around z = 1. Using eq. (2.13) we obtain 



-l-4e 



^14923934 




+ 0{z^) 



, (3.12) 



IBP reduction of the diagrams appearing in the subleading terms gives 




2(l-4e)(3-4e)(l-2e) 



(3-4e)(l-2e) {2e'^ - 2e + 



(3.13) 



(3.14) 



We checked that using these identities we can reproduce correctly the first five terms in 
the soft expansion of Xig. 

The aim of this paper is to compute the leading terms in the soft expansion of the 
triple real emission amplitude for inclusive Higgs production. In order to test our method 
at N3LO, we verified that we can reproduce the correct soft expansion of the phase space 
volume for H + 3 partons. The phase space volume for H + 3 partons in general kinematics 
can be written in the form (see Appendix B) 

1 _5_6,r(l-e)3 



z 



2(47r)5-3e r(6 - 6e) 



2^1(2 - 2e,3 - 3e;6 - 6e;z) 



z 



5— 6€ ^5 



<I>?fe) 



where we defined 



l + (l-e)z + 



me) 



l-e)(3-2e)(4-3^^_,^ 



2(7 - 6e) 
1 r(i - ef 



2(47r)5-36 r(6 - 6e) ' 
Using our method, we obtain the following diagrammatic expansion 



(3.15) 



(3.16) 



^4{z;€) = z 



5-6e 






(3.17) 



All the diagrams in the expansion can be reduced to the soft phase space volume, as 
expected, 






-(1-6) 



(l- 6)(3 - 26)(4 - 36) 
2(7 - 6e) 




(3.18) 
(3.19) 
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To summarize, our method provides a systematic way to perform the threshold expan- 
sion of phase space integrals for the production of a heavy colorless state. Every term in 
the expansion corresponds to a soft integral, as defined in Section 2, which can be reduced 
to a small set of soft master integral using IBP reduction. In the next two sections we 
study some additional properties of soft integrals in general, before applying our method 
to compute the threshold expansion of the triple real emission contribution to inclusive 
Higgs production. 



4. Relations between phase space integrals in different dimensions 

It is well known that loop integrals in different space-time dimensions are related by so- 
called dimensional shift identities [70]. After reduction to master integrals, the dimensional 
shift identities reduce to recurrence relations in the space-time dimension D for the master 
integrals themselves [71-75] . Reverse unitarity rules show that the same dimensional shift 
identities hold for the dual phase-space integrals (see also ref. [76]). 

In this section we present an easy way to derive the dimensional shift identities for 
phase space integrals. We stress that the results of this section are generic and apply to 
phase space integrals in general kinematics. To start, let us consider a phase space integral 
in D dimensions, 

F{D- z^i, . . . , z/„) = j d^N-i{D)f{vi, . . . , z>„) , (4.1) 

where we explicitly indicate the dependence on the space-time dimension D. The integrand 
/ can be written as a product, 

N 

1=1 

where the Pi are polynomials in the rescaled kinematic invariants Sij , raised to some power 
vi. In Appendix A we show that the phase space measure d^N-i{D) for parton-|-parton 
H + [N — 2) partons can be parametrized solely in terms of kinematic invariants, 



d^N-l{D) = Mn-2{D) -,{N-2){D-2)-1 



( \ 

n 

Vv^J,(ij)5^(l,2) / (4.3) 



D-N-l 



X 



N N i-1 ^ 

<5 I 1 - Y^i^u + + -zY.Y.'^A G7^{{sii})Q[GN{{sij})] , 

1=3 i=3 j=3 



where we have defined the normalisation factor 

N 



Mn-2{D) = (-l)^^^^^^2-(^-^)#(27r)(^-^)-(^-2)^ 

1=3 



D 



(4.4) 
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and the Gram determinant Gjv 



GN{{Sij}) = det(siiS2i + SijS2i - Sij)2,<i,j<n ■ (4.5) 

Obviously, the D-dcpcndcnt constants factor out of the integral, and so the actual integral 
depends on D only through the exponent of the Gram determinant. It is then easy to see 
that the phase space measures in shifted dimensions are related by 

d^r,-i{D + 2) = :^(V_2C^+2) _2(Ar-2) d^^_^^D) . (4.6) 

We can thus express a phase space integral in D + 2 dimensions as 

F{D + 2-v^,..., u^) = 1.^^ z^^^-^) f d^N-i{D)GN{{sij})f{ui,...,yn). (4.7) 

For a given set of polynomials {Pi}, the form of the integrand / depends only on the 
exponents {i^i}- One therefore introduces the operators and acting on / as 

(itDi'^l^-- . • • • ,Z^n) = /(J^l,. . . + 1, . . . ,Z^n) , 

_ v4.oj 

ih f ){^l^ • • • , l^i, • • • , Z^n) = /(l^l, . . . , - 1, . . . , f„) . 

If we assume that the {Pi} are linearly independent, we can express the invariants Sij as 
linear combinations of the and /~. This allows us to rewrite the extra power of the 
Gram determinant in eq. (4.6) as a polynomial of degree N in the 

GN{{Sij})f{vi, ...,Un)= Gn{{I-}) fi^l^ • • • ' ^n) • (4.9) 

We thus obtain the following compact formula relating phase space integrals in different 
dimensions, 

F{D + 2; .1, . . . , .„) = ^^^^^^^^.-^(^-^)G^({/r}) HD; -i, . . . , -n) • (4.10) 

Every term in the polynomial can be evaluated according to the action of the /~ operators, 
yielding a superposition of modified integrals in D dimensions. By applying this method 
to a master integral, we can express the master integral in D + 2 dimensions as a linear 
combination of integrals in D dimensions. Using IBP identities, we can reduce the integrals 
in D dimensions to master integrals and thus we find a relation between the master integral 
in D + 2 dimensions and D dimensions. This dimensional recurrence relation can formally 
be written as 

F^{D + 2) = Y,c^j{D)FJ{D), (4.11) 
j 

with coefficients Cij{D) that are determined from the IBP reduction. 
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5. From soft integrals to angular integrals and Mellin-Barnes integrals 



In this section we show that there is a canonical way to derive a MeUin-Barnes (MB) 
representation for the soft integrals that appear in the threshold expansion of phase space 
integrals. The soft integrals we need to consider have the form 

F{e) = J d^ff_,fi{pj};puP2), (5.1) 

where f{{pj};pi,P2) is a ratio of products of multi-particle invariants that is homogeneous 
under a simultaneous rescaling of the final-state momenta, i.e., 

/({Ap,};pi,P2) = f{{Pj};Pi,P2) A'^ , (5.2) 

for some a. Note that this is true only if the (cut) Higgs boson propagator it not raised 
to an additional power. In that case, however, we can always reduce the integrand to a 
homogeneous function using IBP identities. 

Next, we note that there is a subclass of soft integrals that have an additional property: 
they are homogeneous with respect to individual rescalings of the final-state momenta, i.e., 

N 

/({AiPi};pi,P2) = f{{pj};pi,P2) Yi ^7 ' ^^■'^^ 

for some Oj. This subclass of soft integrals is precisely the one where the integrand consists 
of products of power of iiuo-particle invariants, 

m m 

f{{Pj};Pl,P2) = n = Il(^P^. ■PJ.)"'" ' (5-4) 

k=l k=l 

where the index k runs over all the two-particle invariants appearing in /. Every soft 
integral can be converted into an integral of this type, to the price of introducing additional 
MB integrations. Indeed, if we write every multi-particle invariant as a sum of two-particle 
invariants, then we can convert sums into product by using the usual formula 

1 1 /•+^°° dz , , , , , 

_ — r(-z)r(x + z)-—^ , (5.5) 

where the contour separates the poles at z = n from those in z = \ — n, n e N. Without 
loss of generality we can thus assume that our soft integral is homogeneous with respect 
to individual rescalings of the final state momenta. 

If we concentrate on soft integral that satisfy eq. (5.3), it is natural to choose a 
parametrization of the soft phase space that makes the homogeneity explicit. One possible 
parametrization with this property is the so-called 'energies and angles' parametrization, 
where the final-state momenta are parametrized as^ 

PI = ^(1,1,0,...), 

P2 =^(1,-1,0,...), (5.6) 

p.,=^E,(3i, 3<i<N. 



^We work with the rescaled momenta, eq. (2.6). 
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The Ei parametrize the energies of the final-state partons and /3j is the D-velocity in the 
direction pi. In this parametrization the phase space measure for each final-state parton 
takes the form 

d'^pi 6+{pf) = 2-(^-i) e{Ei) -3 dEi dnf-'^ , (5.7) 

where dO,f^~^^ is the measure on the unit sphere parametrizing the solid angle of particle i. 
Furthermore, the on-shell condition for the Higgs boson, coming from the cut propagator 
in eq. (2.16), can be rewritten as 



N 



,P12 ■ {Pl2 - 2p3...Ar 

Thus, the soft phase space measure can be parametrized as 



= S+ {pi2 ■ {pi2 - 2ps,„n)) = 5[l-Y,E^) . (5.8) 



i=3 



TV \ N 



d<I.^_i = (2vr)^-i-(^-2)D 2-(iV-2)(D-i) sIi-^eAH E^-^ dE^ duf''^ . (5.9) 



i=3 / 1=3 



Using this parametrization, we see that every soft integral with an integrand of the form (5.4) 
can be written as 

/m „ m m 

k=l k=l k=l 

(5.10) 

Both the measure and the integrand can be written in a factorized form, and so we can 
integrate out the energies in terms of a generalized Beta function. 

Hence, the only non-trivial integration is a multiple angular integration over the solid angles 
of the final-state partons. Angular integrals can be written in the general form 

In ref. [77] it was shown that such integrals fall into a class of generalized hyper geometric 
functions known as H functions, and an MB representation for the most general angular in- 
tegral of this type was derived. We have thus a general recipe to derive MB representations 
for generic soft integrals. 

Although the previous technique allows us to derive a multifold MB representation for 
every soft integral we need to consider, it can sometimes be useful to insert, if existent, 
explicit closed expressions for the angular integrals. Indeed, for small values of m the 
integrals (5.12) are very simple and can be evaluated in closed form. In the following we 
briefly review some results for angular integrals which will be useful in our case. 

The case m = corresponds to the volume of the solid angle 

^lo-, = J dn^-'^ = • (5-13) 
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Note that this integral is sufficient to compute the soft phase space volume, which corre- 
sponds to putting m = in eq. (5.10): we simply obtain a factor ^d-i for each final-state 
parton. Thus, we obtain 



rfi-e 



9r,-i{e)-[Z7r) l r(2(iV - 2)(1 - e)) ' ^ 

in agreement with the results of Appendix B. 

As we are only interested in massless momenta, /3| = 0, Lorentz invariance implies 
that the angular integral with one propagator must evaluate to a constant. Indeed, we 
have 

~ J iPj ■ P^r ~ r(2-2.-a)- ^^-^^^ 

For angular integrals with two massless propagators one obtains [78] 

^ 22-ai-a2-2e 1-6 r(l - 6 - ai)r(l - € - 0:2) (-5 ^qs 

r(l-e)r(2-2e-ai-a2) 



X 2F1 \ a\,a2\ 1 - e; 1 

To our knowledge, there are no closed formulas for angular integrals with three or more 
massless propagators. We will however need later on the angular integral with three mass- 
less propagators, which admits the MB representation [77], 

^D-i ^ ^ (^Ji ■ P32 5 Ph ' P33 ' Pji ' P33 ) 



r(ai)r(a2)r(a3)r(2 — ai — 0:2 — 03 — 2e) 

f+ioo 

X ' 



j *~ dz^2dzVidZ2^ ^^_^^^^^^_^^^^^^_^_^^^^^^^ ^ ^ Zl3)r(«2 + ^12 + ^23) 



X r(a3 + zi3 + Z23)r(l - ai - 02 - 03 - e - 2:12 - 213 - ^23) 
%i ■ Ph \ ( Ph ' \ ^ ( Ph ' P33 



(5.17) 

with Zij = Zi + Zj , and where the contours separate the poles coming from T functions of 
the form T{. . . — Zij) from those coming from r{. . . + Zij). 

6. Triple real emission phase space integrals in the soft limit 

6.1 Triple real soft master integrals for Higgs production 

After having studied some properties of soft integral in the previous sections, we will use the 
technology developed in the previous sections to compute the threshold expansion of the 
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leading-order cross sections for H plus five partons. More details about the construction of 
the amplitude in this limit will be given in Section 7. Here it suffices to say that we have 
computed the squared amplitude and we have checked that in the limit where we only keep 
the first two terms in the threshold expansion, all the phase space integrals can be reduced 
to linear combinations of the following ten soft master integrals, 





1^ 



2^ 



^2 



{Sl3 + S15)(S23 + S24)S34S35 



^f(e)-^8(e) 



(6.8) 




S15(S14 + Sl5)'S23S34S345 



^f(e)-^9(e) 



(6.9) 
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We have normalized all the integrals to the soft phase space volume for H + 3g defined 
in eq. (3.16). In the remainder of this section we give the dimensional recurrence relations 
satisfied by the master integrals and present the analytic results for each master integral 
as a Laurent expansion in the dimensional regulator e. Technical details about how to 
compute the master integrals analytically will be given in Section 8. 



6.2 Dimensional recurrence relations 

Using the technique described in Section 4, we can derive dimensional recurrence relations 
for all the master integrals defined in the previous section. The knowledge of these recur- 
rence relations provides us with a strong check on our results. In addition, it turns out 
that the master integral J-q{D) is easier to compute in D = 6 — 2e dimensions, where it is 
finite, and the dimensional recurrence relations allow us to relate the six-dimensional and 
four-dimensional results in an easy way. 

The recurrence relation for the soft phase space volume is trivial to obtain from the 
recurrence relation for the T function. 



As we have defined all our master integrals relative to the phase space volume $4 , we can 
simplify their recurrence relations by factoring out the above result. We therefore define 
the ratio 



{D - A){D - 3){D - 2f r(I?-4) 



<l>!iD). (6.11) 



72{D - 1){3D - 5){3D - 4){3D - 2){3D - 1) 6Att^T{D - 1) 






(6.12) 



$f(Z)) 



72{D - 1){3D - 5){3D - 4) (3-D - 2){3D - 1) 



where M was defined in eq. (4.4). We give the results for the remaining master integrals 



relative to TZ. The dimensional recurrence relations for the non-trivial master integrals are 
F2{D + 2)7^ - -3(3^_5)(3^_4^ " 24(3i) - 7)(3D - 5)(3D - 4)-^^^^^ ' ^^"^^^ 



(38 - 2SD + 5L>2) 

MD) , (6.14) 



3(L>-4)(3L>-5) 

(L>- 4)3(1) -3) 



j-4(^ + 2)7^ = - 



18(3D - 10) (3D - 8){3D - 7) (3D - 5)' 
4 (386 - 387D + 1281)2 _ 14d3) 



{D - 4)2(D - 3) 
(£» -4)2(3£> - 14) 



:^4(D) , (6.15) 



J•5(D + 2)7^ 



24(3D - 11) (3D - 8) (3D - 7)" 

(D - 4) (4752 - 9636D + 6706D2 - 1962D3 + 207D4) 
72(D - 3)(D - 1)(3D - 10)(3D - 8)(3D - 5) 
(D-4)2(D-2) 



^5(D) , (6.16) 



Tq{D + 2)7^ = 



96(D-l)(3D-7)(3D-5)' 

(4256 - 6684D + 4224D2 - 1345D3 + 216D^ - 14D5) 



3(D-4)2(D-3)(D-2)2 
) - 10) 
5D-7)" 

(D - 4)3 



(D- 1)(3D-10) 

+ 9(D-2)2(3D-7)-^^^^^ ^^-'^^ 



J-7(D + 2)TZ 



24(D - 2)(3D - 11)(3D - 7) 
4(2D - 7) 



(D-4)(D-3) 

(D-4)'^ 

^72(3D - 11)(3D - 10)(3D - 8)(3D - 7)"^^^^^ ' ^^'"^^^ 



X- I _ 2 (231 - 114D + UD^) 2(D - 4)2(7D - 24) ^ 

^8(^ + 2)7^ - 3p_4)p_3) + 9(D-3)(3D-8)(3D-7)-^'^^^ 

fD-4)'' 

^ ^ -F8(D), (6.19) 



Jb(D + 2)7^ 



72(3D - 11) (3D - 10) (3D - 8) (3D - 7) 
2(3D - 7) (6672 - 7824D + 3460D2 - 684D3 + 51D^) 



3(D-4)2(D-3)2(3D- 10) 

K: 

D 

(D - 4)3 (3D - 14) 



(D-4)(3D-10)(5D-17) ^ , (D - 4) 

12(D-3)2(3D-8) ^^^^6(D-3)(3D-8)^^^ 



96(D - 3)(3D - 13)(3D - 11)(3D - 8) 



MD), (6.20) 



4 (26 - 39D + 16D2 - 2D3) (D - 4)2(3D - 10) ^ 

J-io(D + 2)7^ - (^_4)2p_3) 3(D-3)(3D-8)(3D-7)'^'^^^ 

.^.o(D). (6.21) 



24(3D - 11)(3D - 8)(3D - 7) 
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6.3 Analytic results for the soft master integrals 

In this section wc present the analytical results for the master integrals contributing to 
hadronic Higgs production in the leading and next-to-leading soft approximation. As the 
explicit evaluation of the master integrals is rather long and technical, we delay all details 
about the computation to Section 8 and only summarize the results at this point. The 
first master integral, the soft phase space volume, was already given in eq. (3.16) and will 
not be repeated here. All the remaining master integrals have been evaluated as a Laurent 
series in the dimensional regulator up to terms involving zeta values of weight at most six. 
We have checked that our results agree numerically with the MB integral representation 
for soft integrals derived in Section 5. In addition, the results satisfy the dimensional 
recurrence relations for the master integrals given in the previous section (integrals in the 
shifted dimension have been evaluated numerically using the MB representation) . Finally, 
we make an intriguing observation in our results: if we express all the zeta values up to 
weight six in the basis {C2, Csi C4! C2 Cs) Cs; Cli Ce}) the coefficients in front of the values are 
integers in all cases. We note that this statement is only true in the specific basis of zeta 
values that we chose. The results for the master integrals are listed in the rest of this 
section. 



(Sl3 + S15)S34 

r(6 - 6e)r(l - 2e)2 



3F2(l,l,l-€;2-2€,2-2e;l) 



er(3-6e)r(2-2e)2 
= ^ C2 + 420 Cs - 282 C2 + e (18OO C4 - 1974 Cs + 432 C2) (6.22) 

+ (^5580 C5 + 480 C2 Ca - 8460 (4 + 3024 Qi - 216 C2 



+ (l9260 Ce + 1680 C| - 26226 Cs - 2256 C2 Cs + 12960 (4 - 1512 (3) 
+ 0{e') . 



' S14S23S34 
90 693 1 / \ 1 

^-^ + ^( -6OC2 + I917 +-( - 300 Ca + 462 C2- 2268) 



e2 



- 930 C4 + 2310 C3 - 1278 C2 + 972 e - 2220 (5 - 120 C2 Ca + 7161 C4 

- 6390 C3 + 1512 C2) + ( - 5555 Ce - 300 C| + 17094 (5 + 924 C2 Cs 

- 19809 C4 + 7560 C3 - 648 C2) + O(e^) • 



(6.23) 
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^4 (e) J S13S15S34S45 
3r(6-G()r(l-2e) 
2e4r(l - 6e) 

3r(l - 2e)r(e + 1) 



(1 + 3e)r(l - 3e) 



■ 3 



F2(-3e- l,-2e,-e;-3e,-3e;l) 



+ 4^3(1, 1, 1 - e, -2e; 1 - 2.,! - 2., 2 + 1) 



600 10020 _ 70560 



^ ^480 C3 - 303480) - 3600 C4 + 8016 (3 
- 1007640 - e(^17280C5 - 60120(4 + 56448(3 - 306180o) - 6^(^66000(6 
+ 1920(1 - 288576(5 + 423360(4 - 242784(3 + 918540o) + O(e^) . 

^ J (Sl4 + S15)S23S345 

= (2 - 960 (3 + 684 (2 + € ( - 4620 (4 + 5472 (3 - 1188 (2) 

+ ( - 17160 (5 - 720 (2 (3 + 26334 (4 - 9504 (3 + 648 (2) 

+ ( - 64110 (e - 2880 (| + 97812 (5 + 4104 (2 (3 - 45738 (4 + 5184 (3) 

+ 0{e^). 



1 f d^l 



3 r(6-6e) 



F2(l,l,-2e;l -2e,l-2e;l) 



+ 

e 



^ (^15120 (4 - 59184 (3 + 48600 (2 - 5832o) + 43200 (5 - 207144 (4 
+ 291600(3 - 110160(2 + 23328 + e(lll600(6 - 591840(5 + 1020600(4 
- 660960(3 + 116640(2) + 0{e^) . 



(6.24) 



(6.25) 



^ ^4 (e) y (Sl3 + •S14)(S14 + S15)S23S34 

= ^ - ^ + ^(40(2 + 675) + ^(320(3 - 548(2 - 1530) 

+ ^ (^1500 (4 - 4384 (3 + 2700 (2 + 162o) + 5160 (5 + 320 (2 (3 - 20550 (4 ^^'^^^ 

+ 21600 (3 - 6120 (2 - 648 + e(|l8340 (e + 1280 (| - 70692 (5 - 4384 (2 (3 

+ 101250(4 - 48960(3 + 6480(2) + ©(e^) . 

_ , , 1 r d<^^ 

•^7(e) 



2 e^T{l - 6e) 

■ ^ + ^ + ^ (720 (2 - 24300) + ^ (4320 (3 - 9864 (2 + 55O80) ^g_27) 
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1 r 



^4 (e) J (Sl3 + Sl5)(S23 + 524)534535 

60 822 



e5 



^ (240 C2 - 4050) + ^ (2400 Cs - 3288 C2 + 918o) 
+ ^ (13320 C4 - 32880 C3 + 16200 C2 - 972o) + 51840 Cs + 3360 C2 Cs ^^'^^^ 

- 182484 C4 + 162000 Cs - 36720 C2 + 3888 + e(207600C6 + 11760 C| 

- 710208 Cs - 46032 C2 Cs + 899100 C4 - 367200 Ca + 38880 C2) + C(e 



^9(e) 



J^io(e) 



^4 (e) 7 5i5(si4 + 515)5235345345 
= ^ - ^ + ^( - 120^2 + 2784) + ^( - 120 Cs + 1284C2 + 31968) 
+ ^ (2520 C4 + 1284 Cs - 2088 C2 - 216864) + 15720 Cs + 1920 C2 Cs (6.29) 

- 26964 C4 - 2088 Cs - 23976 C2 + 795744 + e (^82520 Ce + 9600 C| 

- 168204 Cs - 20544 C2 Cs + 43848 C4 - 23976 Cs + 162648 C2 - 244944o) 
+ 0{e^). 

1 f d^f 



^4 (e) J (52s + 524) (524 + 525)534545 

120 2004 14112 1 / \ 
= —^ + —3 ^ + - (240 Cs + 60696) + 1980 C4 - 4008 Cs 

- 201528 + e (6960 Cs + 1680 C2 Cs - 33066 C4 + 28224 Ca + 61236o) ^^"^^^ 
+ (32700 Ce + 6840 C| - 116232 Cs - 28056 C2 Ca + 232848 C4 - 121392 Cs 

- 1837080) + O(e^) . 

7. The threshold approximation to the partonic 2 ^ H + 3 parton cross- 
section 

7.1 General setup 

The production of a Higgs boson in the colhsion of two hadrons hi, h2 is dominated by 
QCD processes. The hadronic cross-section is related to the partonic cross-section by the 
general factorisation formula 

1 

ahi+h2-^H+x = Y^ J dxidx2fi^ixi)fj^{x2)ai+j^HiM'^,xiX2S). (7.1) 

ij 

Here //^(x) are the parton-distribution functions for the parton i inside of the hadron 
h, (Ti^j^H is the partonic cross-section and S is the square of the total centre-of-mass 
energy of the hadronic system. The centre-of-mass energy squared of the partonic system 
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is consequently given by s = XiXjS. The partonic cross-section is expanded perturbativcly 
in the strong coupling constant as- The leading QCD contribution arises in the SM via a 
top-quark loop at 0(a|). 

We consider the Higgs boson to be relatively light compared to the top quark. This 
justifies to work in the limit of an infinite top quark mass and consider Nf light quarks. 
We describe the interaction of the Higgs boson with gluons by introducing the effective 
Lagrangian 

^eff = -\chGI,G'^^''H. (7.2) 

Here G^,^ denotes the gluon field strength tensor and H is the Higgs field. The Wilson 
coefficient ch can be found, e.g. in ref. [79-81]. The next-to-next-to-leading (NNLO) order 
correction to the inclusive Higgs cross-section has been computed in the past by employing 
this effective theory. In this work we present a part of the next term in in the perturbative 
series (N3LO) in the effective theory. 

At every order in perturbative QCD, the cross-section receives contributions from 
various real and virtual radiation processes. We consider only tree-level processes with 
three real-emission partons in the final state. A direct integration over the phase-space of 
the corresponding matrix elements is challenging. As we have discussed in the introduction, 
we pave the way towards the full computation by performing an expansion of the phase- 
space integrals around the kinematic limit where the Higgs boson is produced at threshold, 
s ~ M^. In this limit, all partons emitted in the final state are soft and their momenta 
vanish as z — t- 0. We expand the cross section in the small parameter z defined in (2.5). In 
the following, we shall present the leading and subleading terms in the threshold expansion. 

00 

k=0 

Although we considered in our calculation the Higgs boson H as a final state, we would 
like to stress the universality of our result for the leading term in the soft expansion for 
any other colorless final state produced by gluons in the initial state [82]. The subleading 
term in the soft expansion is no longer universal. 

7.2 Calculation 

To obtain the real-emission cross-section we generate Feynman diagrams using QGRAF [83] 
and compute squared matrix-elements using programs based on GiNaC [84] or FORM [85] 
and our own C-| — h code for color and spin algebra. We perform our calculation in Feynman 
gauge with D gluon polarizations in order to maintain a simple structure for the denom- 
inators of the squared amplitudes. To recover the result for physical gluon polarizations, 
we add matrix-elements with Faddeev-Popov ghosts as external states. We compare our 
results with a set of numeric cross-sections for different phase-space points obtained with 
MadGraph [86] and find perfect agreement. 

In a next step we use rcvcrsc-unitarity and interpret the phase-space integrals as three- 
loop integrals as described in Section 2 and we expand them into a Laurent scries in z. 
The leading and next-to-leading terms in this series are then referred to as the soft and 
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next-to-soft limit of the real emission cross-section. We then derive integration-by-parts 
(IBP) identities for the scalar soft phase-space integrals. We reduce all integrals to a set 
of 10 master integrals by employing the Laporta algorithm [63]. We implemented this 
algorithm in a C-|--|- code that was developed by us specifically for this project, as well 
as the program AIR [87]. The calculation of the remaining 10 integrals is discussed in 
Section 8. 



7.3 Leading soft contribution 



In this section we present the leading soft contributions to the ij^H + X real-emission 
cross-sections, where X where represents the three final-state partons. To obtain a MS- 
renormalised quantity we redefine the strong coupling constant 



as = a§{fiR)e'^^4TrtiRr', (7.4) 



where = r'(l) is the Euler-Mascheroni constant and /jLr is the renormalisation scale. 
In our calculation we consider Nf massless quark flavours. We include a flux factor of ^ 
and average over incoming colours and spins, leading to a factor of 2{i-e)lN'^ -i) ^'^^ gluons 
and 2^ for (anti-) quarks. We abbreviate Cp = ^21^^ ~ -^c- In the case that the 

phase-space integration contains n identical particles we include a symmetry factor of ^. 
For cross-sections with 4 (anti-) quarks we distinguish between the case of only one or two 
quark flavours. The two-flavour case is indicated by explicitly labeling the (anti-) quarks qi 
or Qi. If the (anti-) quarks do not carry any label they are considered all to be of the same 
flavour. We present our results first in terms of the soft master integrals of Section 6. The 
resulting expressions are valid to all orders in the dimensional regulator e. Next, we insert 
the results for the soft master integrals and write the cross-section as a Laurent scries in 
e up to terms containing zeta values of weight six, keeping a factor corresponding to the 
volume of the soft phase-space $4 (e), eq. (3.16), unexpanded. 
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The cross-section containing only gluons is given by 
xi J-i(e)32 ^ 



(3-2e)2e5 (2e2 + e _ l) 

X (336966^^ - 270000e^° + 1341936e^ - 3650136e^ 
+4913370e'^ - 200121 le^ - 3045896e^ + 5040807e'^ 
-33231316^ + 11443306^ - 197535e + 13050) 

64 432e^ + 288e^ - 372€^ - 188e^ + lOe^ + 317e2 - 142e + 15 



+ Me) ^ 

- ^3(e)288 (6e + --7 



- J5(e)32 



e2 (4e3 - 4e2 - 5e + 3) 



e 

126^ + lOe^ + 4e3 - 496^ + 26e - 3 



e2(e + l)(2e-3) 



64 

- -7^7(6)^(6-1) 

64 

-.F8(6)^(e-1) 
^ / x„<,3e2-2e-l 

+ J^9(6)32- 



6e + 1 

2^ 1 1 

3^3!8(iVj~T) 



{A7rasf^f{e)C%CFci, 



218700 2554740 1 / „ ^„ 

— ^ + — — + ^(131220(2 - 9709605) 



+ 4 (782460C3 - 1630854C2 + 14950359) + - (2869830(4 - 9687762(3 



+ +6810588C2 - 8547924) + 8373780(5 + 301320(2 Cs - 35377641(4 
+ 40216932(3 - 11741904(2 + 107996 + e (24995385(6 + 763020(| 
- 103032486(5 - 3541644(2(3 + 145858644(4 - 68849712(3 

+ 7687776(2 - 455984) + 0(e2)|. 

Note that the only colour coefficient of the tree-level matrix element for a Higgs plus 
five gluons is C\Cf (see, e.g., ref. [88]). Obviously this fact is left unchanged by the soft 
limit. 
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In the case that two final-state partons are a quark anti-quark pair we sum over ah 
possible quark flavours and obtain a factor A^^. The corresponding cross-section is then 
given by 



= 8(1 - 6)^(jV2 - 1)2 (4vra5)^$f (6)C^C^4iV^ (7.6) 



xiC 



2 



J'i(e)32 ^ 



(3-2e)2(e- l)e4(e + l)(2e- 1) 
X (288e" - 18288e^° + 129712e^ - 347128e^ + 408738e^ - 107919e^ 
-271140e^ + 359705e'' - 210605e^ + 66680e2 - 10853e + 690) 
^ , , 32 f W8e^ + ISOe^ - 2016^ - Ue^ + 67e^ - 22e + 2 



J^5(e)32 
^i(e)64 



3 V -4e5 8e4 -h e3 - 8e2 3e 

e - 66^ 1 
-2e^ + e + 3. 

12^ - 396e*5 + 982e5 - 1377e^ + 1134e=^ - 527e2 + 122e - 10 



- 3) 



153090 1604043 \ , , , 

— ^3 + ^( - 29I6OC2 + 4903902) 

+ -( - 20412OC3 + 321732C2 - 4833675) - 874800(4 + 2252124(3 - 9IIO88C2 

e ^ ' 

203535 e( - 2711880(5 - 23328OC2C3 + 9651960(4 - 6290136(3 - 492210(2 
1667109) + e2( - 9360360(6 - 816480(| + 29921076(5 2573856(2(3 
- 26589060(4 - 4323186(3 4693212(2 1294731) 

-|- 2CaCf 

321732(2 - 5183163) - 874800(4 2252124(3 - 911088(2 337959 
-f e( - 2711880(5 - 233280(2(3 9651960(4 - 6290136(3 - 492210(2 1651749) 
^ e^{- 9360360(6 - 816480(| 29921076(5 + 2573856(2(3 - 26589060(4 



167670 1743039 \ , n 1 / 

— 4 3 — + ^2 ( - 29160(2 + 5267592) + - ( - 204120(3 



- 4323186(3 + 4693212(2 + 1284491) 



+ 0(e3) . 



In the soft limit only the processes with gluons in the initial state contribute. All other 
contributions containing at least one (anti-) quark in the initial state vanish in this limit. 
This can be understood in terms of soft factorisation: there is no born-level process with 
a massless initial state fermion in Higgs production. 
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7.4 Next-to-leading soft contribution 

In this section we present our results for the next-to-leading term in the soft expansion of 
the real-emission cross sections. Only partonic subprocesses that have at least one gluon 
in the initial state give a non-zero contribution at this order. The results are 
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8. Analytic computation of the master integrals 

In this section we present some details on how to evaluate analytically all the master 
integrals Ti defined in Section 6. The analytic result for the soft phase space volume can 
easily be obtained from the expression for the phase space volume in general kinematics 
and will not be discussed here (see Appendix B for the derivation). 

In general, our strategy is to follow the steps outlined in Section 5: we use the 'energies 
and angles' parametrization to obtain a representation for the integral where the energies 
and angles appear in a factorized form. This may require the introduction of MB inte- 
grations in order to f actor ize sums of invariants in a denominator. We then integrate out 
the energies and angles to obtain a multifold MB representation for each master integral. 
Whenever we are able to do so, we evaluate the remaining MB integrals to all orders in 
e in terms of hypergeometric functions that can easily be expanded into a Laurent series 
e using the HypExp package [89]. In those cases where we did no manage to perform the 
MB integral in closed form for finite values of e, we only compute the Laurent expansion 
of the integral around e = 0, e.g., by resolving singularities in e and summing up harmonic 
sums or by converting the MB integral to a parametric integral which can be computed 
more easily. Details on how to perform these steps for the different master integrals will 
be given in the rest of this section. 

8.1 The master integral T-j, 

In this section we compute the master integral J-3 defined by 



The integrand only involves two-particle invariants, so following the discussion in Section 5 
we can immediately insert the energies and angles parametrization and integrate out the 
energies in terms of V functions. The remaining angular integrations can easily be carried 
out using the formulas given in Section 5. Note that in the present case all the angular inte- 
grals can be performed in closed form, so there is no need to introduce MB representations 
for the angular integrals. We obtain 




(8.1) 



^3(6) = 2' 



,6€-7_3€-3 



r(6 - 6e)r(2 - 2e)r(l - 2e)2 
e2r(2-6e)r(l -e)3 

r(6 - 6e)r(2 - 2e)r(l - 2ef 
e2r(2-6e)r(l-e)3 



J (A-/34)(/32-/33)(Ai-/34) 



= 2' 



,6e-7^3e-3 



(8.2) 




-32- 



The remaining integral can be brought into a more standard form by the change of variable 
cos6'3 = 2y -1, 
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1 (8.3) 
= 2-2^ / dyy-i-^(l-y)-^2i^i(l,l;l-e;2/). 

The integral over y is now easily performed using the recursive definition of the hypergeo- 
metric function 



p+iFp{ai, . . . , Op+i; 6i, . . . ,hp;z) 



T{b,,) 



X I d««^+i-i(l-i)^''-«''+i-%Fp_i(ai,...,ap;6i,...,6p_i;zi). 

Jo 

We immediately get 

^^^'^ = 26^^r(2~-66) ^' 1 - 2^' 1 - 1) • (8-5) 

The 3F2 function can be expanded to the desired order in e using the HypExp package [89], 
and we arrive immediately at the Laurent series of eq. (6.23). 

8.2 The master integral F2 

The integrand of the master integral T2 involves a sum of two-particle invariants in the 
denominator. We replace the sum by a product to the price of introducing a MB integration 
via eq. (5.5), 

d^S f+ioo ^^^^ _ ^ I- ^^s 

'13 ^15 -'34 

The phase space integral is now in the form (5.10), and so we can introduce the energies 
and angles parametrization and integrate out all the energy and the angular variables. This 
results in the following two-fold integral representation for T2-, which is of mixed MB and 
Euler-type, 



/^£l^. p|fir(-.or(. + i) (8.6) 

J (S13 + S15)S34 J-ioo 27rZ J sf3,+ ^ S, .S34 



r(G-G0r(l-2e) , 

1 (8-7) 

xr(l-e + zi) / dyy-' {I- y)-'2Fi {I, zi + 1-1- cy) . 
Jo 

The Euler integral over y could immediately be performed in terms of a 3_F2 function, 
but after that we still need to integrate over the MB parameter zi. We therefore prefer 
not to perform the integration over y, but we rather insert the MB representation for the 
hypergeometric function in the integrand 
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The integral over y evaluates to a Beta function, and we are left with the following two- 
dimensional MB integral 

r(6 - 6f)r(l - 2e) f*'°° dzidz2 , 

r {z2 + 1) r {zi + -22 + 1) r (-e - zi) r (i - e + zi) 

^ r (2 - 2e + Z2) 

The integral over zi is easily performed using Barnes' first lemma, and the remaining one- 
fold MB integral can immediately be recognized as a 3F2 function (see eq. (8.8)). We finally 
obtain the following result for the master integral T2, in agreement with eq. (6.22), 

8.3 The master integral F-j 

The integrand of T'j only contains two-particle invariants, 

$1(6)^^7(6) = / . (8.11) 

We can therefore immediately integrate out the energy and the angular variables. We 
obtain 
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In order to perform the integral over y, we introduce an MB representation for each 2-F1 
function in the integrand and perform the y-integration. This leaves us with the following 
two-fold MB representation for Tt, 
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To proceed, we notice that one of the two integrations evaluates to a 2-^1, which can be 
reduced to V functions using Gauss' identity, 
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Inserting this result into the two-fold MB integral, we immediately see that the remaining 
MB integral evaluates to a 3F2 function, and we get 
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8.4 The master integral 

The integral Fi is defined by 

$f(6)J-4(6)= / . (8.16) 

J S13S15S34S45 

The integrand only contains two-particle invariants and we can immediately integrate out 
the energy and angular variables in the usual way. We obtain a one-fold MB representation, 
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Closing the integration contour to the right and summing up residues at zi = n and 
21 = — 1 — e-|-n, nG , we immediately see that Fa can be expressed as a combination 
of hypergeometric functions, 
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8.5 The master integral 
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The integral F% contains two sums in the denominator, which we can replace by products 
to the price of introducing two MB integrations. 
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We then proceed in the by now familiar way and integrate out the energies and the angles, 
and we arrive at the following two-fold MB representation for J^g, 
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r (-6 + zi - Z2) r (z2 - g) r (-2g -zi + Z2) r (-e - zi + Z2) 

^ r(-e + Z2 + l)r(-2e-zi + Z2 + l) 

Unlike in the previous cases, we were not able to reduce this integral for generic e to 
simple hypergeometric functions. Wc therefore only compute the Laurent expansion of the 
integral. We proceed in the standard way: wc apply the packages MB [90], MBresolve [91] 
and barnesroutines [92] to resolve singularities in e and to expand the resulting integrals 
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under the integration sign and apply Barnes' lemmas in an automated way. The resulting 
MB integrals are at most two-fold, and all of them can easily be done by closing the 
contours to the right and summing up residues in terms of nested harmonic sums defined 
recursively by [93] 

= E ^ S,^n) = ^2^' (8-21) 

k=l k=l 

Note that in the limit n — t- oo harmonic sums immediately reduce to combinations of 
multiple zeta values. The result for reads 
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8.6 The master integral Tio 

The integrand of Tiq involves two sums in the denominator, so we start by introducing 
two MB representations, 
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Integrating over the angles of the particles three and four yields two hypcrgconictric func- 
tions, and we introduce an MB representation for each of them. Performing the integration 
over the last angle, we obtain a four-fold MB representation for Tio- Two integrations can 
immediately be perfumed using Barnes' lemmas, and we obtain 
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After resolving the singularities in e and expanding under the integration sign, all the 
twofold integrals can be reduced to onefold integrals using Barnes' lemmas and their corol- 
laries. The remaining onefold integrals are trivial to compute by closing the contour and 
summing up residues. We find, 
20 2004 141 
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8.7 The master integral 

We start by replacing the sums in the denominator of the integrand of J-^ by three MB 
integrals, 
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We insert the energies and angles parametrization for the final-state particles and introduce 
MB integrations for the angular integrals. Note that the first angular integral necessarily 
involves three massless propagators and can thus not be done in closed form as a 2-^15 so 
we insert the three- fold MB representation (5.17) for it. We then arrive at a representation 
for /"s as a six- fold MB integral convoluted with two angular integrations. It turns out that 
after performing the change of variables ^ Z'i + Z5 and ~^ -^6 + ^1 + ^2) the integrals 
over zi, Z2 and z^ can be done in closed form using Barnes' first lemma. The remaining 
two angular integrations can easily be performed in terms of hypergeometric functions, and 
all but one MB integration can be performed using Barnes' lemmas. We thus arrive at a 
one-fold MB representation for J^^, 

r(6-6.)r(i-3.) 



er(2-66)r(2-2.)r(l-.) 

dzi r {-zi) r {zi + i f r {-zi - 2e) r (zi - e + 1) ^ ' ' 

2m r {zi - 3e + 2) 



One might be tempted to sum up the residues at zi = n and zi = —2e + n, n G N^, for 
finite values of e to obtain an expression for as a combination of two hypergeometric 
functions at 1 valid to all orders in e. The two hypergeometric functions are however 
separately divergent (even for finite values of e) and only their sum is finite. We therefore 
only compute a Laurent series for T^. Resolving singularities in e and summing up residues 
in terms of harmonic sums we obtain, 
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J"5 (e) = —C2- 960 C3 + 684 C2 + e ( - 4620 (4 + 5472 C3 - 1188 C2) 

' ( - 17160 C5 - 720 C2 C3 + 26334 C4 - 9504 C3 + 648 C2) ^g) 
' - 64110 Ce - 2880 C| + 97812 (5 + 4104 C2 Cs - 45738 C4 + 5184 C3) 
+ 0{e') . 



8.8 The master integral Ts 

The master integral J^g is defined by 



$f(e)J-8(e)= / -f^ ^ . (8.29) 

' (S13 + S15)(S23 + S24)S34S35 
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We start by introducing two MB integrations in order to remove the sums in the denomi- 
nator of the integral J-^. The energies are integrated out in terms of F fTinctions, and the 
angular integrations over the particles four and five are performed using cq. (5.16). At this 
stage we have three integrations left to do: the two MB integrations and the integral over 
cos ^3 = 2y — 1 

X r (zi + 1) r {-Z2) r {z2 + 1) r (-e - zi) r (-e - ^2) r (-2e + zi + Z2) ^^'^^^ 

X 2-Fi (1, 2:1 + 1; 1 - e; y) 2^1 (1, 2:2 + 1; 1 - e; 1 - y) . 
In order to proceed, we first apply the identity 

2Fi{a,b;c;x) = {l-x)-^2Fi(^c-a,b;c;-^^ , (8.31) 

and then insert an MB representation for each hypergeometric function in the integrand of 
eq. (8.30). The reason to apply eq. (8.31) before inserting the MB integrations comes from 
the fact that in this way one of the four MB integrations can be performed using Barnes' 
first lemma. We then arrive at the following three-fold MB representation for J^g, 

Gfr(6-G( ) /■+^°° dz2dz3dz4 „ . . „ . . „ . . 

^'^'^=- r{i-errii-6e) L^ ^^-)^r(-2)r(-.3)r(-.4) 

X r (Z3 + 1) r {Z2 - 2e) r {-Z2 - z^) V {z2 + ^4 + l) r (-e - zs) r (^3 - e) (8.32) 
^ r (-26 + Z2- Z3) r (-e - Z4) r {z^ - e) 
^ T (-2e + Z2 + 1) r (-2e - Z3 - Z4) 

In the rest of this section we show how we can compute a Laurent expansion for this 
integral. We proceed in the standard way and resolve singularities in e. At the end of 
this procedure, we have a collection of MB integrals of dimensionality at most three with 
integration contours that are straight vertical lines. These integrals can then be safely 
expanded in e under the integration sign. In the following we discuss the computation of 
the two and three-fold integrals. 

Three- fold MB integrals. There is one three-fold integral contributing to J^s, 

L ' (~^2) r {Z2) r {-zif r {za) r {-Z2 - za) r {-z^ - z^) 

^'"^ , (8.33) 

^ r (1 - Z3 - zj) r {z2 + Z4 + 1) r (zg + -^4) r {z2 + zz + z^) 

r(z2 + i)r(z3) 

where we omit all T function prefactors and where the integration contours are straight 
vertical lines defined by 

Re(z3) = 0.28 and Re(z4) = 0.97 and Re(z5) = -0.36 . (8.34) 
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We start by closing the Z3 contour to the right and take residues, 

r+'^dz2du r {-Z2) r (^2) r (-^2 - z^) r {-z^f r (^4) r {z2 + ^4 + 1) 

J-ioo (27ri)2 T{z2 + l) 

T{n + Z2) T{n + Z2) T (n + Z2) 

-1p [n + Z2) rV' [n - Z4) 




(8.35) 



riT (n — z^) nT (n — Z4) F (n 



The first two sums can be performed in terms of F functions and their derivatives. We 
ihustrate this on the first term (the second term is similar), 

^ T(ji + Z2) . d ^ T{n + Z2 + r]) 

/ — W/ \r [1^ + Z2) = hm — > — — — 

■^nF n — Z4) v-^o Of] ^—i nF (n — za) 

n=l ^ ' n=l ^ ' 

r;^o dr] y F (-Z4) J (8 36) 



r(^2) 



F(-Z4) 



■0 (2:2 + 1) V' (-22 - -2:4) - — (--^2 - -2:4) - Z2ii {Z2 + 1) V' (--24) 

^2 



+ —V' (-2:4) - -0' {-Z2 - Z4) 
Z2 

In this way, the first two terms can effectively be reduced to the computation of two-fold 
integrals, and we will therefore not discuss them any further in this section. 

The third term can also be summed up in closed form. However, unlike the first two 
terms, the sum cannot be expressed in terms of F functions and their derivatives alone, but 
the sum evaluates to a 4F3 function. We therefore arrive at the following single three-fold 
MB integral, 

r+'^ dZ2dz4 F {-Z2) F {Z2) F {-Z2 - Zi) F {-Z^f F (Z4) 

y-.oo (2vri)2 F(l-Z4) (8.37) 

X F (Z2 + Z4 + 1) 4i^3 (1, 1, 1, Z2 + 1; 2, 2, 1 - Z4; 1) . 
Although it looks like we have managed to reduce the three-fold integral to a two-fold 
integral, it is still secretly three-fold, except that we have 'hidden' one integration inside the 
4F3 function. The advantage of this representation is that we can change the representation 
for the 4F3 function in the integrand. More precisely, we perform the change of variables 
24 ^ — Z4 — Z2 and chose the contours to be straight vertical lines given by 

Re(z2) = - and Re(z4) = - . (8.38) 
3 5 

We then insert an Euler integral representation for the 4^3 function as an integral over a 
3F2. The 3F2 function turns out to be reducible, 

3F2(l,l,l;2,2;t) = (8.39) 

and so we finally arrive at 



7^8,3 =-r^ ^ f dte^-^ (1 - t)--^ Li2(t) 
J-ioo K^^V Jo 

^ r {-Z2) F {z2) F (1 - Z4) r (-Z2 - Z4) F {z2 + Z4) 

X ', 

r Z2 + 1 



(8.40) 
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Next, we would like to exchange the MB and the Euler integration and sum up the residues 
of the poles of the T functions in the integrand. However, we are only allowed to do so if 
the Euler integration does not produce any new poles whose residues need to be taken into 
account. It is easy to see that in our case the Euler integral converges whenever Re{z2) 
and Re(^;4) are positive. We can thus close both contours to the right and exchange the 
Euler and MB integrations and then sum up the residues coming from the F functions. 

We start by taking residues in Z2 and then the residues in Z4. There are several cases 
to be considered separately: 

1. The poles at Z2 = n2 € give rise to poles in Z4 at 714 G N^. We hence obtain 
double sums of the form 

n2 + nA S^^{n2) St^in^) St^{n2 + n^) 



"-2 /t4 -^hV^A) Ot:i{ll'2^ 'lH) Q _ /g 

^ ^ n2 J ni' ni' {n2 + n^^^ ^ ' ' ^ ' ' 



n2,n4=l 

Sums of this type can be performed using XSummer [94, 95] , and give rise to compli- 
cated multiple poly logarithms whose arguments are rational functions of t and {\ — t). 
Using symbols (see Appendix D), all these complicated multiple polylogarithms can 
be reduced to harmonic polylogarithms with indices and 1 in t. 

2. Taking the residues at the poles at Z2 = — -2:4+712, n2 G gives rise to the expression 

,\tu2{t)± ^(1-0-^ (8.42) 

J-ioo 27r« (rz2 - Z4) V (712 - 24) 

Next we want to close the z^ contour and sum up the corresponding residues. From 
the previous discussion, we are forced to close the contour to the right, and the 
summand obviously only has poles at ^4 = 714 G inside the integration contour. 
There is however a subtlety, and we cannot just sum up the residues. Indeed, it is 
easy to see that 

• for 714 < ^2, there are double poles. Shifting the summation variable n2 — > 
n2 + ?T-4, these residues give rise to double sums similar to eq. (8.41) and can again 
be performed using XSummer. We obtain complicated multiple polylogarithms 
with rational functions of t as argument. Using symbols, they can again be 
simplified to harmonic polylogarithms with indices and 1 in t. 

• for 7i4 = n2, there is a simple pole. This gives rise to a simple sum which is 
trivial to perform in terms of harmonic polylogarithms in t. 

• for 714 > 7Z2, there is a simple pole. Shifting the summation variable 714 — >■ 712+714 
this gives rise to a single double sum S{\ — t, 1 — 1/t), with 

"2 ,"-4=0 

This sum is not of the type (8.41), and wc therefore need a different way to sum 
up the series. This procedure will be discussed in the rest of this section. 
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One way to sum the series S{x, y) is to recognize that it is related to an integral over an 
Appell function. More precisely we can write 



S{x,y) = V 



n2!n4! x^^+i y"4+i 



xy ^ ^ (n2 + ^4 + 1)! ^2 + 1 n4 + 1 

n2,n4=0 ^ ' 



= _± r,, V -Miw^ q x; (8.44) 

7o h n ("'2 + ^4 + 1)! "2! ^^4! 

n2,n4=U 

= -— / rdxF3(l,l,l,l;2;e,x). 
The particular Appell F3 function we obtain turns out to be reducible, 

F3(a, 7 - a, /3, 7 - 7; x) = (1 - x)"+^~^ 2i^i(a, 7; C + X - Cx) , (8.45) 
and so we obtain a simple two-fold integral representation for 5, 

S{x,y) =-— fd^ r dx2Fi{l,l;2;^ + x-^x) 

xy Jo Jo f . 

1 r^. /■^^^, iog(i-g)+iog(i-x) _ ^- ^ 

a^y 7o io ? + X - ?x 

This integral is trivial to perform in terms of multiple polylogarithms with rational func- 
tions in x and y as arguments (see Appendix D). It then follows that S'(l — t,l — 1/t) can 
be expressed in terms of multiple polylogarithms with rational functions in t as arguments. 
Using symbols, we arrive at the following simple expression. 



5 { 1 - t, 1 - = -^^^-^ [ - 4Li3(t) + 2Li2(t) log t 

+ Jlog3t + ^logt + 4C3 
6 3 



^.47) 



Putting everything together, we arrive at a representation for J-8.3 as a one-fold integral 
over harmonic polylogarithms. This integral is trivial to perform, and we obtain 

Two-fold MB integrals. There is only one two-fold MB integral contributing to J-g 
that cannot be reduced to simpler integrals by Barnes' lemmas and their corollaries. This 
integral reads 

_ , , , -'^ dz3dzi r i-zsf r (z3) r (^3 + 1) r (-^4)' r (^4) r (^4 + 1) 

•>^8,2(e) 



-lOO 

X 



(27ri)2 2er (-Z3 - Z4) (8_49) 

3etp {-Z3) + eip {Z3) - 2eV' (--23 - ^4) + etp (-2:4) + eip (24) - 1 , 



where the integration contours are straight vertical lines defined by 

Re(z3) = -0.64 and Re(z4) = -0.22 . (8.50) 
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The integral could in principle be done by closing contours to the left and summing up 
residues. This leads to double sums of the form (8.41) - thousands of them due to the 
presence of multiple poles - but without any parametric dependence. We therefore took a 
different route, which we present in the following. 

We start by noting that the integrand of eq. (8.49) agrees, up to higher order terms of 
0{e), with the function 

r (^3 + 1) r {-Z4f T {z4 + 1) r (-e - ^3)' r (^3 - e) r (-e - ^4) r {z^ - e) 

2er(-2e - zs - z,) " ^^'^^^ 

It would however be wrong to conclude that then necessarily the integrals, seen as a Laurent 
expansion in e, are also equal to the same accuracy, because the Laurent expansion of 
cq. (8.51) around e = might require to shift the integration contours to avoid pinch 
singularities in the limit e — )■ 0. It is however easy to see that, for the contours given in 
eq. (8.50), no new pinch singularity is created for e ^ in eq. (8.51). We therefore consider 
from now on e to be infinitesimal but finite: it is large enough to separate poles at, e.g., 
Z2 = from Z2 = e, but small enough to ensure that no poles change their nature, i.e., poles 
that were left (right) of the contour (8.50) in eq. (8.49) remain left (right) of the contour 
in eq. (8.51). If we chose e in this way, we conclude that 

T8,2{e) = ^8,2ie) + 0{e), (8.52) 

where ^8,2 is given by eq. (8.51) integrated over the straight vertical lines defined by 
eq. (8.50). Our aim will be to find an e expansion for J^8,2- 

We perform the change of variables {z3,Z4) (—1 — z^,—Z4), and we close the z^ 
contour to the right and take residues. There are two towers of poles we need to take into 
account: 

Z3 = ns G N and Z3 = -1- e + ns, ns G . (8.53) 
Two comments arc in order: 

1. At this stage our assumption that e be infinitesimal but finite is vital, because other- 
wise the two towers of poles would glue together and thus give rise to double poles. 

2. The second tower of poles runs over the set { — e, — e + 1, . . .}. For technical reasons 
that will become clear below, it is easier to explicitly take into account the residue 
at ^3 = — 1 — e, and to compensate for this by adding it back, 

^8,2(e) =7^8,2(e) + F8,2(e), (8.54) 

where ^8,2 is the integral obtained from J^8,2 by deforming the Z3 contour such that 
the pole at Z3 = — 1 — e is now to the right of the contour. The residue at 2:3 = —1 — e 
can be computed in closed form and gives rise to 

V (A- r(i-6)Mi + ^)r(i-26)3 . + 2-n 

7^8,2(6) i6,5(i + ,)r(l-36) 1' 1 - 1 - ^ + 2, 1) 

r(l-2e)^r(l-e)2r(l + e)^ _ 2>T{l-2efT{l-efT{l + ef 
^ 8e6r(l-4e) 16e6r(l-3e) 
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Next we compute Fs^2 by closing the contour to the right and summing up residues. The 
resulting sums can easily be performed in terms of hypergeometric functions, and we get 



i^8,2(e) = - / 



+^°° r (1 - Z4) r {z4f r (-e - Z4) r {z4 - e) 



2iTi 2eT {zA - 3e) V (-2e + Zi + \) 
r(-e - l)r(l - efV {Z4 - 3e) 3F2 (1 - e, 1 - e, 1 - e; e + 2, -2e + ^4 + 1; 1) 

+ r(-2e)^r(e + l)r (-2e + Z4 + I) 3F2 (-2e, -2e, -2e; -e, 24 - 3e; 1) 

(8.56) 

In order to proceed, we insert an Euler integral representation for each of the functions. 
It is then easy to see that the Euler integrals are convergent for Re(z4) > and e infinites- 
imal but finite, and so we can exchange the Euler and MB integrations provided that we 
close the integration contour to the right. The important point is that the 2F1 functions 
appearing inside the Euler integrals are independent of Z4, and so they can be pulled out 
of the MB integral. Summing up the residues in 2:4, we then arrive at the following integral 
representation for ^8,2- 

FsM^) = ^^^y^ j\u-'{l-t)-' (8.57) 
r(i - 2e)2r(l + e) 



4e(l + e)t 



2Fi(-2e, -2e; -e; t) 2^1(6 + 1, e + 1; e + 2; 1 - t) 



, ,r(l-e)^r(l + e) 



. r(i-6) 3 
(i + e)^ 



,\ , 2i^i(l - e, 1 - e; e + 2; i) 2Fi(e + 1, e + 1; e + 2; 1 - t) 



The terms involving a single 2-^1 function in the integrand immediately evaluate to 3F2 
functions, which can be expanded in e using HypExp. In addition, the fourth term can 
be done in closed form as follows: We insert an MB representation for each 2-^"! and 
perform the Euler integration as a Beta function. One of the two remaining MB integrals 
evaluates to a 2-^1 evaluated at 1, which reduces to T fTinctions through Gauss' identity. 
The remaining one-fold MB integral then immediately evaluates to a 4F3 function, which 
can be expanded in e using HypExp. We were not able to find a closed form for the last 
remaining Euler integral. We therefore insert an Euler integration for each 2-P^i) and we 
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obtain the expression, 

r(i-6mi + 6) 

2e3r(l + 2e) ^^'^ ' ^^'^^^ 

where I denotes the integral 

X(e) = / dtdudv{l-t)-'{l-uf'u-'v'{l-tuy-\l-{l-t)v)-'-^ . (8.59) 
Jo 

It is easy to see that X is finite as e — )■ 0, and so we can expand in e under the integration sign 
and perform the integration over t, u and v recursively using the techniques of Appendix D. 
This is a trivial exercise that leads to 

= 2C, - |le + (25& - ^6) + (-6C| - + . (8.00) 

We have thus obtained the e expansion of ^8,2, and thus of Fs,2- We find, 

^8,2(e) = (^y Ca - ^Cs) ^ - 187EC5 + 4C| + ^T^^'Ca - ^ + 0(e) , (8.61) 
where = r'(l) denotes the Euler-Mascheroni constant. 

The result for Ts- We have now computed all the two and three-fold integrals con- 
tributing to Ts- The remaining one-fold integrals are trivial to compute and we obtain 

Tsie) =-^ + ^ + ^(240C2 - 4050) + ^(2400(3 - 3288(2 + 918o) 

+ - f 13320 C4 - 32880 G + 16200 C2 - 972o) + 51840 (5 + 3360 C2 G 

e ^ ^ (8.62) 

- 182484 C4 + 162000 Cs - 36720(2 + 3888 + 6^207600(6 + 11760(1 

- 710208 Cs - 46032 (2 Cs + 899100 C4 - 367200 C3 + 38880 C2) + C(e^) . 
8.9 The master integral Xg 

In this section we describe the computation of the most complicated master integral, the 
integral 

$f (6) Me) = [ , ^ "^^f . (8.63) 

We start in the usual way and derive an MB representation for J-g by inserting the energies 
and angles parametrisation and inserting MB integrations for the angular integrals. After 
applying Barnes' lemmas and their corollaries several times, we arrive at the following MB 
representation for J^g, 

J^9(e)=J^9,i(e)+J^9,2(e), (8.64) 

where 

3(1 + 6e)r(6 - 6e)r(l - 4e) f+'°° dzidz2dz3dz. . , 

•^^'^-2e(i + 4e)r(i-6e)r(i-e)3y_ i2m)^ r(-.,)r(-.3) 
X r (1 - zi - 2:2) r (1 - zi - Z3) r {zi + Z2 + Z3) r {-zi - Z2- z^- z^) 

(8.65) 

xV{zx + Z2 + Zi + l)T {zi + Z2. + Zi + l)T {-zi -Z2- 2e) 
r (zi - e - 1) r (-Z1 - ^4 - e - 1) r (z4 - e + 1) 



r (1 - Z2) r (1 - Z3) r {-zi -Z2- 4e) r {z^ - 2e + 1) ' 



-44- 



_ 3(l + 6e)r(6-6e)r(l-4e) [+'°° dzidz2dz3dz^ . 



2e(l + 46)r(l - 6e)r(l - 6)3 J_^^ {2^^ 

X r (-Z3) r - Z3 + 1) r {Zl +Z2 + Zs) r {-Zl - Z2- Zs- Z^) 

X r (zi + Z2 + + 1) r (zi + Z3 + Z4 + 1) r (-zi - 2:2 - 2e) 

r (zi - e - 1) r (-^1 - Z4 - e - 1) r (2:4 - e + 1) 
^ r (1 - 22) r (1 - za) r (-zi - 4e) r (^4 - 2e + l) " 

We can resolve singularities for J^g^i and J-9,2 and expand in e up to and including 0{e). 
The result is a collection of high dimensional MB integrals which, after closing the contour 
and taking residues, result in multi-fold harmonic sums. While we were able to perform all 
the harmonic sums in terms of zeta values for all MB integrals up to 0{e^), at 0{e) new 
polygamma functions appear in the integrand which make the combinatorics of the sums 
rather intricate. We therefore chose a different method to evaluate the integral J^g, which 
we describe in the rest of this section. 

We start by noting that the J-g is finite in Z? = 6 dimensions. This can easily be 
checked by replacing e by e — 1 in the MB representations (8.65) and (8.66) and resolving 
singularities. Our goal is to find a parametric integral representation for J-g in D = 6 — 2e 
dimensions and to expand under the integration and perform the parametric integrations 
recursively. The result in D = 6 — 2e can then be related to the (divergent) result in 
D = A — 2e using the dimensional recurrence relation for /"g of Section 6. 

It is easy to derive a parametric representation for using the technique described 
in Appendix C. We find 

. X r(12 - 6e)r(3 - 3e)r(l - e) r ^ .J 

T^{D = 6 - 26) = r(5-66)r(2-6)4 F'^^'^ +^9,2(6)J , (8.67) 

with 



-1 +l-2e 



X9i(e) = - / dtidt2 / dxidx2dx3tl~^'' {1 + tiY 
Jo Jo 

X x^' (1 - xif~^' xl~^' (1 - X2y' x^ (1 + ^2^3)^"''^' (1 + ^2:^2x3)' (8.68) 

X (titlxiX2X3 + 1^X2X3 + tit2XiX2 + tit2X3 + 12X2X3 + t2 + h + l)^^ ^ , 



Xg,2(e) = / dhdt2 / dxidx2dx3tl-^^l + tiy-Ul-^^ 
Jo Jo 

X x]-' (1 - xif-^' x\-^^ (1 - X2Y' X3 ^ (1 + t2X3f~'^' (1 + 12x2x3)' (8.69) 

X {t\t^X\X2X3 + t^X\X2X3 + t2X\ + t\t2X\X2 + t\t2X3 + t2X\X2X3 + tl + Xi)' 



v3e-3 



Several comments are in order about the parametric integrals we just defined. First, one 
can easily check that both Xg^i and Xg^2 are individually finite as 6 — >■ 0. Second, at first 
glance our goal to integrate out the integration variables one-by-one seems rather hopeless 
due to the appearance of the huge polynomial factor. However, as we will see shortly, 

there is a sufficient condition that allows one to test whether a parametric integral can 
be performed in terms of multiple polylogaritlims, and this criterion is fulfilled for the 
integrands of Xg^i and Xg^2- We very briefly summarize this criterion in the following, and 
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we refer to ref. [96] or to Appendix D for more details. In order to understand the criterion, 
it is important to first understand multiple polylogarithms and their integration. 

Multiple polylogarithms are generalizations of the ordinary logarithm and the classical 
polylogarithms, 



dt dt 

— and Lin{z) = J yLi„_i(t), 



(8.70) 



with Lii{z) = —ln{l — z). Multiple polylogarithms are defined recursively via the iterated 
integral [97,98] 

P dt 

G{ai,. . . ,a„;z) = / G(a2, . . . , On; t) , (8.71) 

Jo t - ai 

with G{z) = 1 and where ai,z G C. In the special case where all the a^'s are zero, we 
define, using the obvious vector notation a„ = (a, . . . , a). 



G{6n; z) = - In'' z. (8.72) 
n! 

In the special case where G {—1,0, 1}, multiple polylogarithms are related to harmonic 
polylogarithms, 

H{ai, ...,an;z) = {-If G{ai, . . . , a„; z) , (8.73) 

where p denotes the number of elements in (ai, . . . , a^) equal to +1. 

Suppose now that we are given some integral over some rational function. We would 
like to be able to perform the integrations one-by-one in terms of multiple polylogarithms. 

Unfortunately, such a procedure is not possible for a generic integral. In ref. [96] a suffi- 
cient condition was derived that allows one to test whether such a recursive integration is 
possible. For a detailed description of the criterion we refer to ref. [96], or to Appendix D, 
where (a version of) the criterion is reviewed. In a nutshell, as can easily be seen from 
the definition of multiple polylogarithms, a sufficient condition to perform the parametric 
integrations recursively using multiple polylogarithms is that at each integration step we 
can find an integration variable in which all the denominators appearing in the integrand 
are linear. More formally, let S^'^^ be the set of all polynomial factor in the integrand 
of Xg^j. Next, we define the sets S^Jy y where x,y,... are integration variables, as the 
sets of irreducible non-monomial polynomial factors in the integrand after integration over 
{x, y,---) in this order. For more details how to construct these sets we refer to ref. [96] and 
to Appendix D. Here it suffices to say that the sets S^J^ ^ can be constructed solely from 
the knowledge of the sets S**^*-* . A sufficient criterion for a recursive integration in terms of 
multiple polylogarithms is that there is an ordering of the integration variables such that 

(i) 

the sets S^J ^ ^ corresponding to this ordering contain at least one integration variable in 
which all the polynomials of the set are linear. Miraculously, this condition is satisfied for 
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both Ig^i and Iq^2- We have 



= {1 + t^, 1 _ xi, 1 - X2, 1 + t2X3, 1 + t2X2X3, 

tit\xiX2X^ + t\x2X^ + tit2XiX2 + tit2X^ + it2a;2a;3 + ^2 + ii + 1} , 

5(2) = {1 + i-^^ 1 _ a;-^^ 1 _ a;2, 1 + t2aJ3! 1 + ^23^23^3, 

tit\xiX2X^ + t|xiX2X3 + ^23^1 + tit2XiX2 + tit2X'i + t2a;ia::2a;3 + + 

^l^^^^ = {1 - Xl, 1 - X2, 1 + ^23^1X2, 1 + t2a;3, 1 + t2X2X'i, 
- t2X^X2 + t2XiX^X2 + X1X2 - X3X2 + X3 - 1} , 

(2) 

Sl^^^^ = {1 - Xi, 1 - a;2, 1 + t2XiX2, 1 + ^22^3, 1 + i2X2X3, 
t2Xi - t2X2Xi + t2X2X3Xi - t2X-i + Xi - 1} , 
1 - X2, 1 + t2X2, 1 - X3, 1 - 
t2X2X3 + X2X3 - X3 + 1} , 



•^(tl.xi) = {1 - a;2, 1 + ^2X2, 1 - X3, 1 + i2X3, 1 + ^2X2X3, 



(2') 

S[t[,xi) = {1 - 3:2, 1 + ^2^2, 1 - X3, 1 + t2Xi, 1 + t2X2Xz, 



- t2X2 + t2X2X3 + ^2 + 1} , 



^'(tl,xi,X2) = {1 ~ ^3, 1 + hXs, 2t2X3 -t2 + X3} , 
= {2 + t2 - 3:3, 1 - X3, 1 + t2Xs} , 
S^L,..,.,)={^ + t2,l+2t2}, 



8.74) 
8.75) 

Xl} 

8.76) 

8.77) 

8.78) 

8.79) 

8.80) 
8.81) 
8.82) 
8.83) 



We see that if we perform the integration in the order (ti, xi, X2, X3, ^2) then at each step 
ah the polynomials are linear in the next integration variable. The actual integration can 
be carried out in an algorithmic way. The procedure is however rather lengthy, so we do 
not discuss it here in detail, but we refer to Appendix D for a detailed description of the 
integration algorithm. The result is 

J='g{D = 6- 2e) = 1663200C3 - 55440071^ + 3326400 
+ 120c ^ + 13097r'' - 244203(3 + 286l7r^ + 135294^ 

- 2e^ j - 25779600C5 - 9702007r^C3 + 8386577r'' - 8149392(3 - 2017567r^ 

(8.84) 



- 31378284j + ^9605757r^ + 180873000C| - 1978358850(5 

- 336120757r2C3 + 56663280(3 + 324050l7r^ + 68361307r^ + 810381510^ 
+ O(e^) . 

Using the dimensional recurrence relations for J-g derived in Section 6 we then finally find 
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the value of J^g in D = 4 — 2e dimensions, 

160 1712 1 / \ 1 / \ 

J-9(e) =^ ^ + ^ - 120 C2 + 2784 j + ^ - 120 Ca + 1284 C2 + 31968J 

+ ^ (2520 C4 + 1284 Cs - 2088 C2 - 216864^ + 15720 (5 + 1920 C2 Cs 

- 26964 C4 - 2088 (3 - 23976 C2 + 795744 + e ^82520 Ce + 9600 C| ^^"^^^ 

- 168204 Cs - 20544 C2 Cs + 43848 C4 - 23976 C3 + 162648 C2 - 244944o) 
+ 0(e'). 

9. Conclusions and outlook 

In this article, we presented a method for the expansion of phase-space integrals in kine- 
matic parameters. We applied our method to perform a threshold expansion of triple-real 
radiation partonic cross-sections for Higgs production at hadron colliders. Our method 
reduces the phase-space integrals in the coefficients of the expansion to a small number of 
soft master integrals. These master integrals emerge not only in Higgs production processes 
but in all processes which are topologically equivalent, such as DrcU-Yan production. We 
have computed the master integrals using a combination of techniques for Mellin-Barnes 
integration and the integration of real parametric integrals. We presented explicit expres- 
sions for the first two coefficients in the threshold expansion of all triple-real radiation 
partonic cross-sections which are necessary for the inclusive Higgs cross-section at N3LO, 
both in terms of soft master integrals and as a Laurent series expansion in the dimensional 
regulator e. 

We can apply our method to compute further coefficients in the threshold expansion. 
This requires the reduction of higher rank cut loop integrals in their powers of propagators 
and irreducible numerators. It also requires the evaluation of some new master integrals 
originating from diagrams which did not contribute to the first two terms of the threshold 
expansion. We believe that it is feasible to compute in this way a sufficient number of 
terms in the threshold expansion for phenomenology purposes. A second possibility is to 
apply the reverse-unitarity method to compute the same cross-sections without resorting 
to a threshold expansion. The work that we have presented here will be particularly 
important for solving the unexpanded master integrals. With reverse unitarity, one can 
derive differential equations for the master integrals in arbitrary kinematics. For their 
solution, we can use the soft limits presented here as boundary conditions. Subleading 
terms in the expansion can serve as a check of the solutions or as an aid to guess them. 

A complete computation of the N3LO coefficient requires also virtual contributions in- 
tegrated over the phase space of up to two real partons in the final state. The corresponding 
one-loop, two-loop and three-loop amplitudes are already known in the literature as a Lau- 
rent expansion in e [99-109]. Integrating them over phase-space requires in addition their 
universal single-real and double-real infrared limits at higher orders in e. A threshold ex- 
pansion of mixed real and virtual corrections requires an extension of our method. While 
the threshold expansion for real emission matrix elements can be performed by uniformly 
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rescaling all final-state parton momenta, pi — >• zpi, and expanding in the small parameter 
z, for real- virtual matrix-elements the scaling of the loop momentum cannot be determined 
without further consideration. To determine the correct expansion at threshold we consider 
the strategy of regions [110,111]. We applied this method to all 2-loop real- virtual master 
integrals contributing to Higgs production at NNLO and could reproduce the soft limits 
of the integrals. In the following we will illustrate this technique on the example of the 
integral 



2> 



{P3 +P2)^ 



k^{k + pif{k - pzf{k + P12Y ■ 



(9.1) 



(9.2) 



We introduce the light-cone coordinates, e.g. for g'^, 

and denotes the remaining {D — 2) Euclidian transverse components. Furthermore we 
choose our reference frame such that 



Pi, 



P2, 



71' 

0, P2,- = 



0, Pi,± = 0, 



1 



(9.3) 



V2' 



PK 



0. 



We make the soft region manifest by redefining the loop momentum and the final state 
parton momentum in terms of the soft-expansion parameter z 

k-^zk-, k±^V^k±, p^^zps- (9.4) 

After expanding in z we find the leading soft integral to be 



-l-3e 



1 



d'^k- 



PS+ J 'A;2(P + fe+)(A;2 + fc+p3_)(A;+ + l)' 

All external kinematic scales are factored out of the integrand and the integral can be 
easily computed using Schwinger parameters. This yields 

T(l-e)2r(6 + l) 



-3 2 



Z 



in agreement with the result in ref. [67]. In full kinematics other regions, which are sup- 
pressed in the soft limit, contribute to the integral. 

Finally, the computation of the collinear counterterms for the partonic cross-section has 
been performed in refs. [112,113]. Alternatively, this task can easily be achieved numerically 
following the procedure of ref. [36] and using the results for the partonic cross-sections at 
NNLO through 0{e) of refs. [67,68]. 

We believe that the prospects for a complete computation of the N3LO coefficient 
for Higgs production and other processes are excellent, and we are looking forward to 
completing this task in future publications. 
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A. Parametrisation of the 2 ^ H + {N — 2) phase space 

In this appendix we derive the parametrisation of the 2 — >■ i7 + (A^ — 2) phase space in 
terms of Lorentz invariant quantities. The usual Z?-dimensional phase space measure is 
given in eq. (2.3) and we rescale the momenta as defined in eq. (2.6) and find for s = 1, 

d^N-i{pH,P3, ■ ■ ■ ,Pn;z;D) 

(A.1) 

for two incoming partons with momenta p\ and p2 producing a Higgs boson with momentum 
Ph and n additional partons with momenta Pi. We parametrise the final-state momenta 

as Pi = {Pifi-,Pi,z-,Pi,±), so that Pi^±_ is the {D — 2)-dimensional transverse component of 
the momentum. Then we can rewrite the integration over the momenta of the final-state 
partons as, 

dPpi5+ (Pi) = dpifldpi^zd^~'^Pi,±@ {pi,o) S {pIo - p^^ - pj^) , (A.2) 
Next we choose a frame such that 

pi = ^(l,l,0,...,0) and p2 = ^ (1,-1,0,..., 0) . (A.3) 



In this frame we can write the f-channel invariants {i > 3) as 

Sli = {Pl - Pif = -'^PlPi = - {Pi,0 - Pi,z) , 



S2i = {P2 - Pif = -2p2K = - {Pi,Q + Pi, 



(A.4) 



or equivalenty 



such that 



Pifi = ~ {sii + S2i) and pi^^ = ~ {s2i - su) , (A.5) 



pI± = suS2i ■ (A.6) 
The s-channel invariants {i, j > 3) can be written as 

Si,j = isiiS2j + sijS2i) - 2pi^±_ ■ pj^± . (A. 7) 

This yields 

d^Pi6+ [pf) = ^d^~'^pi^±dsuds2i@ (-sii) & {-S2i) 5 {siiS2i - pl±) ■ (A.8) 
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After eliminating the integration over the Higgs momentum using momentum conservation 
the measure becomes 



i=3 — (A.9) 



It is well known that rotationally invariant integrals of the form 

i 

can be rewritten by making the integration over the scalar products manifest: 

^ ^ 2^ / n ^D+i-i n d{Pi-Pj) [Gram(pi, . . . ,pn)] ^ Q [Gram(pi, . . . ,pn)] fiPrPj) ■ 

i=l 0=1 

(A.11) 

Here Gram(pi, . . . ,pjv) denotes the Gram determinant of the vectors pi, which is defined 
as 

Gram(pi, . . . ,pn) = det(pi ■Pj)i<i,j<N- (^-12) 

We observe that the measure as well as the matrix element are invariant under rotations in 
the transverse space. Therefore, we can use this relation to rewrite the D — 2 dimensional 
transverse integration 

N N i 

n d''~''pi,± = 2'-"^ n ^^+1-* n ^(^^.^ • ^'^•.^) 



i=3 1=3 j=3 

X [Gram(p3,_L, ■ ■ ■ ,Pn,±)]^^ © [Gram(p3,_L, ■ ■ ■ ,Pn,±)] 



(A.13) 



Using the definition (A. 7) of the s-channel invariants, we can rewrite the integration over 
the scalar products as 

n n d{Pi,± ■ Pj,±) = (-2)^^^^"^ ({[ dpl^ ( n ^s,, ] . (a.i4) 

i=3j=3 \i=3 J \3<i<j<N J 

The integration over the norm of the transverse components can be performed using the 
on-shell (5-functions, so that we find 

N 

(jV-2)(iV-3) c o„r (iV-2)(JV-3) 



Ud^PiS^iP^) = (-l)^^^^^2«-3^--^ 



i=3 

/ N 



D+l-i 



( \ 

n 

l<i,j<N 



(A.15) 
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where we have written the Gram determinant as 

GN{{sij}) = det{siiS2j + SijS2i - Sij)3<ij<N ■ 

Thus we finally arrive at the phase space measure 



(A.16) 



n ^^i- 



\<i,j<N 



N N i-1 

i=3 j=3 



1=3 



with 



N 



i=3 



as in eq. (4.4) and 



D 



D 



(A.17) 



(A.18) 



(A.19) 



B. Derivation of the phase space volume 

In this appendix we derive an expression for the 2 — > + (iV — 2) phase space volume that 
is valid for all orders in z and e. We start by recalling the general phase space factorisation: 

If we assume that we know the phase space volume at N'^LO, i.e. for 2 — > + A;, 

$fe+i(m2,s) = j d^k+iim'^,s), (B.2) 

this relation allows us to rewrite the phase space volume at N'^+'LO as 

^k+l+i{m\s) = J J ^^d^i+i{fi\s)d^k+i{m\fx^). (B.3) 

We are specifically interested in rewriting the phase space volume at N^+^LO as a convo- 
lution of the N'^LO phase space with a two particle phase space. Specialising to Z = 1, we 
can use this formula to inductively derive the phase space volume for arbitrary orders. We 
have, 

f f '^ du? 

$fc+2(m2,s)=/ / — d$2(//^s)c^$fe+l(m^M2) 

J Jrn? 

= I ^d$2(/i^s)$fe+l(m^/x2). 



(B.4) 



-52- 



Prom eq. (A. 17) we obtain an explicit parametrisation of the two particle phase space and 
we find 

X <^(1 - Xi^k+3 - X2,k+3) ~) ^k+i{m'^,fJ''^)- 

We can perform the integral over and X2,k+3 in terms of beta functions and make 

the transformation /j? = sx, = sz to obtain 



^k+2{z, s) = (47r)-^+^5^-^ ;;; J' dx (l - xf-'^^^k^^^zs, xs). (B.G) 



r(i-e 

r(2 - 2e 

In the following we use eq. (B.6) to prove inductively the following result: 



1 rn — f\"- 

^ 2^ ^ r(2n(l-e)) ^ ' (B.7) 

X 2F1 ((n - 1)(1 - e),n(l - e),2n(l - e); 1 - z) , 

First, eq. (B.7) correctly describes the phase space volume for n = 1 and n = 2. In order 
to derive $n+2 iteratively from $^+1 we use eq. (B.7) and find 

$ ^oiz S) = l('47r)l-2("+l)+("+l)^ r(l - e)"+^ n-{n+l)e 

'^n+2[z,s) ^[^w) r(2n(l - e))r(2 - 2e) 



Z \ 2n— 1— 2ne 



X dx{l - a;)i-2^x"-i-"^ (1 - £ 
X 2F1 ((n - 1)(1 - e), n(l - e), 2n(l - e); 1 - ^) . 
To solve the integral, we make the transformation a; = 1 — (1 — z)y and find 

$„+2(z, S)=C [' dy y'-^'il - y)2"-l-2ne^i _ _ ,)y)-n+ne 

Jo 

X2F1 f(n-l)(l-e),n(l-e),2n(l-e);(l-z) ^"^ 



(B.8) 



1 - (1 - z)y 
where we have factored out 



(B.9) 



f _ ^(A \l-2(n+l)+(n+l)e ^(1 e)"+^ ri-(n+l)ef-, _ x2(n+l)-l-2(n+l)e (y, -, 

^-2^^''> r(2n(l-e))r(2-2e)'^ ^ • ^ ^' 

Next, we introduce a Mellin-Barnes representation for and exchange the MB integration 
with the parametric integration. This allows us to perform the integration over y in terms 
of another 2-^1. Then wc find 

* ,)-c r(2n(i - 6))r(2 - 26) f+^^ dj , r((n - i)(i - 6) + 

^^'^ ' ^" r((n-l)(l-e))r(n(l-e))y__ 2^^^' ' ^ r(2(n + 1)(1 - e) + 
X r(n(l - e) + O2F1 (n(l - e) + 2 - 2e, 2(n + 1)(1 - e) + 1 - z) . 

(B.ll) 
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Introduce another MB integral for the we arrive at 

r((n - 1)(1 - 6) + Or(n(l - 6) + ^ + x)r(2 - 2e + x) 
r(2(n + l)(l-e) + e + x) 

In the next step we perform the change of variables C ^ C + X) such that 

r((n-l)(l-e))r(n(l-e))7_ioo (27r«)^ 
r((n - 1)(1 - e) + C - x)r(n(l - e) + 0^2 - 2e + Q 
r(2(n + l)(l-e)+e) 

The integral over x can now be performed using Barnes' first lemma, and we find 



r((n - 1)(1 - e))r(n(l - e)) 27ri 
r(n(l-e) + Or((n + l)(l-e)+0 ^ ' ^ 

r(2(n + l)(l-6)+0 

The integral over ^ is just the MB representation of a 2F1 so that we finally find 

$ ^ofz - lr47r)l-2("+l)+("+l)S"-("+^)"ri - ~)2(n+l)-l-2(n+l)6 ^(1 - e)"+^ 

*n+2i^,sj- 2^ s [L z) r(2(n + l)(l-e)) 

X 2F1 (n(l - e), (n + 1)(1 - e), 2(n + 1)(1 - e); 1 - . 

(B.15) 

We have therefore inductively shown the validity of (B.7) for all n. 

C. Prom Mellin-Bcirnes integrals to parametric integrals 

In this appendix we describe how one can derive an Euler-type integral from a Mellin- 
Barnes integral with a balanced integrand. Roughly speaking, a MB integral is said to be 
balanced if for each integration variable Zi the number of T functions of the form r(. . . — Zi) 
is equal to the number of F functions of the form T(. . . + Zi). More precisely, the integral 



/. 



^llTiak,+Zir'=il[Tibk,-Zif>'^, ak,,/3k,eZ, (C.l) 



ki = l k2 = l 



is said to be balanced if Yl^^ = Y1^2 ' assume in the following that the contours 
are straight vertical lines such that the real parts of the arguments of all the F functions 
are positive^. In that case we can always derive an Euler-type integral representation for 



Note that in dimensional regularisation we might need to require e to be finite for such a contour to 
exist. 
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the MB integral. We start by noting that if an integral is balanced, then we can always 
express its integrand clS Si product of Beta functions, 

B{x,y) = / + (C.2) 

Jo 

The integral (C.2) is convergent whenever Re(x),Re(y) > 0. It is easy to convince oneself 
that this condition is satisfied whenever the real parts of all arguments of the F functions 
are positive in the original MB integral. We can therefore replace each Beta function by 
its integral representation (C.2) in the integrand of the MB integral and, because all the 
integral are convergent, we can exchange the MB integrations and the integrations coming 
from the Beta functions. This leaves us with an integral of the form 

/ n ^o(t) Reity / n v^Rii^r- ' (c-3) 

Jo / J-ioo ^TTZ 

where t = (ti, . . . ,tjv) and the Rk are ratios of products of the ti and 1 + tj. Next, we 
would like to perform the MB integrations. This can be done using the formula 



zo+ioo J 

a^ = S{l-a), a > . (C.4) 



27r? 



Indeed, parametrizing the contour as z = zq + it, we obtain 

rzo+ioo 

' ZQ—ioo 

Equation (C.3) can thus be written in the form 

N \ M 



/ J^a' = a'o / — e^*^'^« = a^o5(lna) = (5(l-a). (C.5) 

Jzn-ioo 27rz 2tt 



r [tl^A ^o(t) Reity n 5 (1 - Ri{t)) . (C.6) 

•^0 \n=l / m=l 



We can solve the S constraints, and the result is the desired parametric integral. 

In the following we discuss two very simple examples that illustrate the above proce- 
dure. We start by discussing Barnes' first lemma, i.e., we consider the integral 



+'°° dz 

— r{a + z)T{b + z)r{c-z)r{d-z). (C.7) 

ZTTl 



We assume that the integration contour and a, b, c and d are such that the real parts of 
all r functions are positive. We rewrite the integrand in terms of Beta functions, 

r+i°° dz 

I = r{a + c)r{b + d) B{a + z,c-z)B{b + z,d-z) (C.8) 

/+ioo J foo 
^ / dtidt2tl+'-\l+ti)-''-Hl+'-\l+t2)-^''^ 
-ioo 27r« Jq 

POO 

= r(a + c) T{b + d) / dti dt2 (1 + 4"^ (1 + hy''^'^ S{1 - iiis) , 

^0 
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where the last step follows from eq. (C.4). Solving the (5-function constraint leads to a 
one-fold integral that can immediately be recognized as a Beta function, and we recover 
the usual form of Barnes' first lemma, 

POO 

I =T{a + c) r{b + d) j dti tl+'^-^ (1 + 

_ r(a + c) r(a + d) r{b + c) r{b + d) 
~ r{a + b + c + d) 

The second example we are going to discuss is Gauss' hypergeometric function. We 
consider the integral 

We again assume that all conditions for convergence are satisfied. Rewriting the integrand 
in terms of Beta functions, we obtain 



T(b) /'+*~ dz 
J=—±±- B{b + z,-z)B{a + z,c-a)x' 

r(c - a) J_i^ 2m 



_V{b) 

r(c 



^ r dt, dt2 tr' (1 + ii)-^ 4"' (1 + 12)-' s f 1 - . 

-a) Jo \ 1 + 



(C.ll) 



Solving the (5-function constraint with respect to t2 and performing the change of variables 
ti — >■ ~ 0) immediately arrive at the usual integral representation for the 2-P1 
function, 

J = J^^^ , r dh n-' (1 + h)'-' (i+h+x h)-' 

L[c — a) Jq 

r(c - a) Jq 

= -p . \ {a, b; c; -x) . 

r(c) 

D. Symbolic integration of parametric integrals 

In this section we describe an algorithmic approach to compute certain classes of para- 
metric integrals. We stress that this approach is not genuinely new, but has already been 
successfully applied, in some variant or another, to the computation of Feynman integrals, 
e.g., [96,114-118]. More precisely, consider an integral of the form 

where and 6^. are integers and the Pfc({xj}; {yj}) are polynomials with integer coefficients, 
which wc assume irreducible over Z, i.e., they cannot be factorized into a product of non- 
constant polynomials of lower degree. We assume that the integration range is [0, 00] for 
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each integration variable Xj. While many statements remain true for generic integration 
boundaries, in certain cases the algorithm we are going to describe breaks down, due to the 
appearance of boundary contributions. We can however always map a generic integration 
region Xi G [a, b] to yi G [0, oo] by the change of variable 




(D.2) 



We furthermore assume that the integral is convergent for e = 0, and so we can expand in e 
under the integration sign. Note that in the simplest divergent cases where the singularities 
in the integrand factorize we can reduce the problem to a convergent integral by subtracting 
the divergencies. 

Our goal is to compute the coefficients of the Taylor expansion of I{{yj}; e) by integrat- 
ing out the integration variables recursively one-by-one in terms of multiple polylogarithms, 
eq. (8.71). Obviously, not every integral can be performed in this way. In the rest of this 
appendix we describe a sufficient condition, first obtained in ref . [96] , to be satisfied by the 
integrand of I{{yj};e) so that it can be integrated in terms of multiple polylogarithms. 
We then review an algorithm for integrating these classes of integrals, and illustrate the 
procedure explicitly on a simple example. 

D.l Denominator reduction 

In this section we review (a variant of) the sufficient condition of ref. [96] to determine 
whether an integral can be performed recursively in terms of multiple polylogarithms. The 

main idea is that we need to determine an ordering of the integration variables such that 
at each step during the integration all the denominators are linear in the next integration 
variable. We start by defining the set S of all the polynomials that are not monomials and 
that appear inside the integrand of eq. (D.l), 



To start the integration, we have to assume that there is one integration variable, say Xa 
such that all the element of S are linear in Xa- In that case we may write 



where Q%{{xi}; {yj}) and Rf.{{xi}; {yj}) are polynomials that arc independent of Xa- Note 
that, after expansion in e, the integrand may also contain logarithms of the Pk{{xi}; {yj}), 
which can be rewritten in terms of multiple polylogarithms. 




(D.3) 



Pk{{^^h {%}) = QUi^ih {%}) + Rl 



{{xi}\{y3]) 



(D.4) 



\ogPk = \ogiQlxa + Rl) = \ogRl + \og[l + ^x, 




(D.5) 




where for clarity we suppressed the arguments of the polynomials. Furthermore, we can 
use the shuffle algebra of multiple polylogarithms to replace every product of multiple 



polylogarithms by a sum, 



G{ai, . . . ,ani; z) G{ani+1, ■ ■ ■ ,ani+n2'i — G'(<^(t(1)5 ■ ■ ■ )0(^(ni+n2); (D 6) 

a-eS(ni,n2) 

where S(ni,n2) denotes the set of all shuffles of ni + n2 elements, i.e., the subset of the 
symmetric group Sni+n2 defined by 

S(ni,n2) 

_i _i _i _i (D-7) 

= {cr G Sni+nJcT ^(l)<...<cr ^(rii) and cr ^(rii + 1) < . . . < cr ^(ni+n2)}. 

Thus we can assume without loss generality that the integration over Xa takes the form 

f°° dxa 

Jo (Qtxa + . . . + R^J-^r. ^(«; • (D.8) 

Partial fractioning the factors in the denominator, e.g., 

(D.9) 



(Qpa + RtKQfxa + Rf) QtRf - QfRt K^a + Rt/Qt Xa + Rf/Qf 

we obtain a sum of integrals that can be reduced to the recursive definition of multiple 
polylogarithms, eq. (8.71). We can thus easily compute a primitive with respect to Xa, and 
then take the limits Xa ^ and Xa ^ oo. We will address the issue of limits in the next 
subsection. 

We would like to iterate this procedure and integrate over the next variable. This is 
however only possible if inside the new integrand we can still find an integration variable in 
which all polynomials are linear. The polynomials appearing inside the integrand are 
and R^, which have been introduced through eq. (D.5) and (D.9), as well as the combina- 
tions Q'^Rf — QfRk introduced by partial fractioning. This last polynomial however need 
not necessarily be linear, even if the and i?^ are. In order to proceed, it is therefore 
mandatory that all the Q^Rf — QfR^ factor into polynomials that are linear in a certain 
variable. 

In ref. [96] a criterion was given that allows one to determine a priori whether the 
above procedure terminates, i.e., whether there is an ordering of the integration variables 

such that all the denominators stay linear at each integration step. We start by defining 
the set 5'(a,^) as the set of irreducible factors that appear inside the polynomials Q'^, i?^ and 
Q%Rf — QfRk- Then, if we can find an integration variable Xf, such that all the elements 
of S(^xa) linear in Xb, we can restart the above procedure and integrate over Xf,. If we 
iterate this procedure and are able to construct a sequence of sets of polynomials 

S{Xa)^ S(xa,Xb)^S(xa,Xb,Xc)^ ■ ' ' (D.IO) 

such that in each set all the polynomials are linear in at least one integration variable, then 
we have found an ordering of the integration variables such that wc can recursively integrate 
out all the integration variables in terms of multiple polylogarithms. We stress that this 
condition is sufficient, but not necessary: even if we fail to find a suitable sequence (D.IO), 
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the integral might still be expressible in terms of multiple polylogarithms (e.g., after a 
suitable change of variables). In addition, note that for the last integration step it is not 
necessary for the polynomials to be linear: we can then factor the polynomials into linear 
factors, whose roots involve algebraic expressions of the parameters {yj}- 

D.2 Symbolic integration 

We have now a criterion to determine if a given integral of the type (D.l) can be inte- 
grated recursively in terms of multiple polylogarithms, and we have already explained how 
to perform the first integration step. We assume from now on that we have found a se- 
quence (D.IO) and the associated ordering of the integration variables, and that we have 
performed the integration over Xa as described in the previous subsection. 

In order to continue and perform the next integration, we need to address two issues 

1. We need to take the limits — >■ and Xa — >■ oo of the primitive with respect to Xa- 

2. If we want to integrate next over Xh, we have to face the problem that Xb might be 
'hidden' inside the arguments of multiple polylogarithms of the form 

«(...,-|....). (D.U, 

i.e., Xb appears inside the polynomials and R'j^. In order to compute a primitive 
with respect to Xb, we need to rewrite all multiple polylogarithms of the form (D.ll) 
as G{a; Xb), where a is independent of Xb- 

The limit Xa — )■ can easily by taken by using the fact that 

lim G(ai;xa) = 0, ai ^ 6 . (D.12) 

While eq. (D.12) is sufficient to compute the value at Xq = in many circumstances, it can 
happen that the primitive has spurious poles and / or logarithmic singularities at Xa = 0. 
In the presence of poles we need to expand the multiple polylogarithms around Xa = 0. 

This task can easily be achieved using the series representation of multiple polylogarithms. 
Indeed, multiple polylogarithms can equally well be represented as multiple nested sums, 

^^mk,...,mi{Zk, ■ ■ ■ , Zi) = ^ ■ ■ ■ •^^W ~ -^m2,...,TOfc(?^ ~ 1;^2, • • • , Zk) , 

ni>...>nfc>0 ^ ^ n=l 

(D.13) 

where -^m2,.--,mfc('^ ~ '^y Z2, ■ ■ ■ , z^) denote the Z-sums defined in ref. [94]. The integral 
representation and the series representation are related by 

G(0^_^,ai,...,0^_^,afc;x) = (-l)'=Li^,,...,^, (^^,...,^^ , ^ . (D.U) 

mi — I TUk — l 

If Ofe = 0, we can use the shufHe algebra to extract pure logarithms in Xa, e.g., if 6 7^ 0, 

G{b, 0; Xa) = G(0; x„)G(6; x„) - G(0, b; x„) . (D.15) 
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We can thus easily obtain the expansion around Xa = in terms of Z-sums to arbitrary- 
order. 

Next, we discuss how to take the hmit Xa oo. If we define Xa = ^/xa, then the 
problem can be reduced to taking the limit — > 0, 



We would like to use the inversion relations for the multiple polylogarithms in the right- 
hand side to express them as linear combinations of multiple polylogarithms of the form 
G {. . . ;Xa), for which the limit Xa ^ is trivial. Unfortunately, the corresponding func- 
tional equations are often unknown. Note at this point that the problem of finding the 
inversion relations is formally equivalent to the second problem we have not yet addressed, 
namely how to bring the primitive with respect to Xa into a form where Xf, only enters 
through multiple polylogarithms of the form G{a;xb): in both cases we are looking for 
functional equations that bring a certain multiple polylogarithm into a 'canonical form' 
where a certain variable only appears as the explicit upper integration boundary of the 
multiple polylogarithms. In the rest of this section we argue that, under certain conditions 
which we will later show to be equivalent to the criterion derived in the previous subsection, 
plus the hypothesis that the integration ranges are [0, oo], we can always bring multiple 
polylogarithms into the 'canonical form' 



for some variable x such that is independent of x and the coefficients Cj involve only 
multiple polylogarithms that are independent of x. Note that this result is similar to the 
result obtained in ref. [96]. In the following we give a constructive algorithm that allows 
us to derive the canonical form (D.17). 

In order to achieve a rewriting of our multiple polylogarithms in canonical form, we 
need to derive the corresponding functional equations. The natural language to discuss 
functional equations among multiple polylogarithms are symbols [119-123] and the Hopf 
algebra of multiple polylogarithms [98]. We start by giving a concise review of symbols. 

One possible way to define the symbol of a multiple polylogarithm is to consider its 
total differential [98], 




(D.16) 




(D.17) 



dG{an-i 



n-l 



. . . , a^-i, Oj+i, . . . , ai; a„) d\n 



( ^1^^) , (D.18) 
\ai- Qi-i ) 



and to define the symbol recursively by [122] 



n-\ 




(D.19) 



S{G(an-\ 



,---,ai;an)) = ^5(G(i 



. . . ,ai_i,aj+i, . . . ,ai;an)) 
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As an example, the symbols of the classical polylogarithms and the ordinary logarithms 
are given by 

S{Un{z)) = -(1 - z) (8) z (8) ■ ■ ■ (8) z and S In" = . . ^.®z^ . (D.20) 

(n— 1) times n times 

In addition the symbol satisfies the following identities, 

...(g) (a - ft)®... = ...®a® ... + ...®h® ... , 

(±1)(8)... =0, (D.21) 

S (Cia; x) G{h- y)) = S{G{d; x)) in S{G{h- y)) , 



where III denotes the shuffle product on tensors. 

We can make the following observation about the symbol of a multiple polylogarithm: 
if the Oj are independent of x, then the symbol of G(ai, . . . ,a„;x) contains exactly one 
term hat contains x in all its entries, and this term can be chosen to be of the form 

5(G(ai, . . . , a„; x)) = (a„ - . . . (ai - x) + . . . . (D.22) 

In order to proof this statement, we focus on the term in the total differential of G(ai, . . . , a„; x) 
proportional to dx, 

dx 

dG{ai,...,an;x) =G{a2,...,an;x) h... 

x-ai (D.23) 

= G(a2, . . . ,a„;a;)dlog(ai - x) + . . . , 

where the last step follows from 

, X dx — dai „ ,s 

dlog(ai-a;) = -, (D.24) 

X — ai 

and where the dots indicate terms in the total differential that are independent of dx. The 
statement then follows recursively. Note that this statement is independent of whether the 
ai are zero or not. 

Assume now that we are given an integrable symbol T (which will correspond later to 
the symbol of the multiple polylogarithm we want to bring into canonical form) which is 
of uniform weight w and has rational coefficients. If T does not satisfy this last condition, 
we deal separately with the contributions of different weight and / or different rational 
function prefactor. Let us suppose that the entries are drawn from a set S. Without loss 
of generality, we may assume S to consist of irreducible polynomials in some variables Xi, 
1 < i < n, and for simplicity we assume for now that the polynomials are linear in all 
the Xi (we will see in the next subsection that the correct restriction is that S satisfies the 
reduction criterion of the previous subsection) . Furthermore, we assume that we have fixed 
an ordering on the variables, which we will take in the following as {xi, . . . ,Xn). For each 
variable Xj, we define a linear map (^a;. which acts on elementary tensors s by 

^^^(^) ^ I C (-^, . . . , -f^; X,) , if s = (a,,x, + 6^) ® . . . ® (aix, + h) , ^^^ ^5) 
, otherwise . 
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Morally speaking, the map (j)xi assigns to T the combination of multiple polylogarithms 
which will give the same terms that have Xi in all entries through eq. (D.22). 

Using the maps ^3;., we can now formulate an algorithm that assigns to T a multiple 
polylogarithm in the canonical form associated to the ordering of the variables (xi, . . . , x„). 
We start by defining a new symbol by subtracting off the contribution from (px^ (T) , 

Ti=T-S[(l)x,{T)]. (D.26) 

By construction, each term in the symbol Ti has at least one entry that is independent of 
xi. Next, concentrate on the terms of the form 

XI (Hi,...,in,bii^{ai2Xi + bi^)^...^{ai^xi + biJ), (D.27) 

where by hypothesis the 6^^ are independent of xi. As we have subtracted of the contri- 
bution from (pxiiT), these terms cannot come from a multiple polylogarithm of the form 
G{. ■ ■ ;xi) of weight w, but it can only arise from the product 

logbi, G ( . . . , - ) ^ (pxAK) <Px, {{ai^xi + bi^)(^.. .(^{ai^xi + bij)) . (D.28) 

It is easy to convince oneself that the difference 

T2 = Ti- X Cii,...,i„5 (px2ibh)4>xi(^{ai2Xi + bi^)0 ...0{ai^xi + bij)j (D.29) 

contains only terms for which at most {w — 2) entries depend on xi. We can now go on an 
recursively subtract contributions with different multiplicities of Xi . Assume for example 
that we have subtracted all contributions where Xi appears in more than {w — r) entries, 
and that the resulting symbol is T^. We can then concentrate on the terms 

XI bi^® ...®bi^® {^ir+iXi + bi^^^) ® {ai^xi + 6^^)) . (D.30) 

It is easy to convince oneself that in the difference 
Tr+i = Tr— 

X ^i,-,iw ^ ^X2 (bil <8) . . . (8) (l>xi (^{air+iXi + 6i,+i) (g) . . . (g) (a^^Xi + kj)^ 

{ii,...,iw) 

(D.31) 

xi appear in at most {w — r — 1) entries. We continue this procedure until we reach Tyj, 
which is independent of xi, and we restart the algorithm with and (I)x2- We then repeat 
this procedure until we have exhausted all the integration variables, and the algorithm 
stops. The result of this algorithm is, by construction, a function of the form 

X Cii,...,j„G(ai„;x„)...G(aji;xi), (D.32) 

(ii, •••,»«) 
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whose symbol is T and such that the are sequences of rational functions that are in the 
variables Xk, k > ik, i.e., the sought-for canonical form for T. 

Using this algorithm we can easily take the limit — >^ oo, by applying it to the right- 
hand side of eq. (D.16) with the ordering {xa,Xb, . . .). We can also use it to rewrite the 
integral over Xa such that we can easily integrate over the next integration variable Xf,. 
However, the result of the algorithm will at this stage still deliver the wrong answer, as 
we have so far only constructed a function in canonical form whose symbol matches the 
symbol of a given function. To illustrate this, assume we run the algorithm on the function 
in the right-hand side of eq. (D.16) with the ordering {xa,Xf,, . . .). The result is a function 
Q{xa-, xjj,. . .) in canonical form such that 



It would however be wrong to conclude that G [a; — j and Q{xa, Xf,,. . .) are equal, because 
the symbol maps to all terms proportional to multiple zeta values. In ref. [124] an algo- 
rithm was described that allows to reconstruct these zeta- valued terms using the full Hopf 
algebra structure of multiple polylogarithms, augmented by some ideas by Brown [125]. 
In the following we only give a very brief account on how to reconstruct the zeta-valued 
terms, and we refer to ref. [124] for a detailed description of the algorithm. 

In the following we assume for simplicity that the function G ^a; — G{xa, Xb, . . .) is 
reaF, so we do not need to worry about imaginary parts proportional to itt. Next we act 
with A2,i,...,i on the difference, where A2,i,...,i is the component of the iterated coproduct 
where the first component has weight two and all other components have weight one. 
Without loss of generality we can write 



where the a^j, are irreducible polynomials and Ai^ is a combination of multiple polyloga- 
rithms of weight two. Without loss of generality we can assume that we have collected all 
term that have the same 'tail' logOjj . . • <8) logai^_^, i.e., we can assume 

logai2®...(8)logai„_, / log O . . . (g)log aj^_i if (i2, . . . , / (j2, . . . • (D.35) 

As we know that the symbol of the function vanishes, eq. (D.33), we necessarily conclude 
that is proportional to ^2, i.e., = ki^C,2^ for some rational number fcj^. This rational 
number can easily be determined by evaluating Ai^ numerically at a single point using any 
of the standard libraries to evaluate multiple polylogarithms [84,126-131], and running for 

There is no obstacle to consider complex-valued functions. Imaginary parts can be extracted in exactly 
the same way using Ai,...,i. 




(D.33) 





(D.34) 
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example the PSLQ algorithm [132]. Equation (D.34) then takes the form 



A. 



2,1,.. .,1 



G [a;—] - g{xa,Xb, ■ ■ ■) 

Xa , 



E 



_i hi C2 <8) log a 



12 



log ai^_^ , 



(D.36) 



(il,...,iu,_l) 



Next we drop (2, i-e-, we only keep the tail of each elementary tensor. If we also drop the 
log signs, we obtain a symbol associated with the terms proportional to (2, 



A 



2,1,., 



X, 



(n,...,iTO-i) 



_l (X- 



'»2 



(D.37) 



Running the algorithm described at the beginning of this section on the symbol in eq. (D.37) 
we obtain a function Q2{xa-, x^, . . .) of weight u; — 2 in canonical form such that 



A2,l, 



G [a;^ ] - g(xa,Xb,.. •) - C2 ^2(a;a, a;^, 



0. 



(D.38) 



We have in this way determined all the contributions proportional to (2, and the result 
is by construction in canonical form. We then repeat exactly the same exercise by acting 

with A3J 1 to determine the terms proportional to and we continue in this way until 

we have exhausted all possibilities and the algorithms stops. As a result we obtain the 
expression of G (ji; J-^ in canonical form at function level. 

This terminates the algorithm to bring multiple polylogarithms into the canonical form 
corresponding to a certain ordering of the variables. This operation is sufficient to take the 
limit — 7- 00, and to bring the integral into a form where the next integration can easily 
be performed. We have however not yet shown that this algorithm will always terminate 
on the class of integrals we consider. In particular, when formulating the algorithm, we 
assumed that all entries that appear in the symbol are linear in all variables. In the next 
subsection we show that this condition can be relaxed, and that the algorithm always 
terminates if the integrand satisfies the reduction criterion of the previous subsection and 
if the integration range is [0, 00]. 

D.3 Symbolic integration and denominator reduction 

In this subsection we show that the algorithm we just described always terminates for inte- 
grals of the type (D.l) such that the set S of polynomials satisfies the reduction criterion. 
Prom now on we assume that we have found an ordering of the integration variables, which 
we take as (xi, . . . , x„) and we can find a sequence 



(D.39) 



such that all the elements of <S'(a,^ .^fe) linear in x^^i. We start by showing that under 
these hypotheses and after having integrated out (xi, . . . , x^.i), 

1. the symbol of the primitive with respect to x^ has all its entries drawn from the set 



^ (xi,...Xk-i) ~ {Xki • • • 1 Xn} U S'j-j.j^ 



(a:i,...a;fe_i)- 
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2. the symbol of the function after integration over x^, i.e., after taking the limits 
Xfc — )■ 0, oo of the primitive, has all its entries drawn from S'(a;j,...a;j, ) = {xk+i, ■ ■ ■ , x^jU 

S{xi,...x^)- 

We start by proving the first statement, the second then immediately follows by taking the 
appropriate limits. As in the following we constantly switch between polynomials of the 
sets we just defined, and polynomials as entries of symbols, we introduce the notation [P] 
to refer to P as an entry of a symbol. Note that we then have the identities 



[PQ] 
[P/Q] 
[±i] 

etc. 



[P] + [Q] , 
[P] - [Q] , 

0, 



(D.40) 



We proceed by iteration in the number k of variables we have already integrated out. 
We start by analyzing what happens after the first integration. It is obvious from eqs. (D.5) 
and (D.9) that after the first integration the primitive only involves multiple polylogarithms 
G{ai , . . . , ; xi ) with 

aie{0,-Rl/Qi} . (DAI) 

It is then easy to see (e.g., from the polygon approach to the symbol of ref. [123]) that the 
symbol can only have the following entries: 



[xk],[Rl], 

[-Rl/Ql]=[Rl]-[Ql], 

[Rl + Qlxi]-[Ri] 



1 + 



1 + 



xiQl' 
Pi - 

QjPi 
PlQ] 



[Ok P] - Pk Qi] 



= m - [pi] 
[pi] - [Qh ■ 



(D.42) 



The polynomials that appear inside the symbol are precisely those that appear in S and 
'^{xi)j which finishes the proof of the first statement for the first integration. 

Next consider taking the limits xi — >• and xi oo. By definition, Rj. and Qj. are 
independent of the limit, so we only need to consider the limits of [xi] and [Pk]. [xi] will 
give rise to logarithmic singularities in the limit, and these terms must cancel if the integral 
is convergent. For [Pk] we have. 



limfPfc] 
Um [Pk] 



[Pk] ' 

Ihn {[x,Rl + Ql] 

xi—>-0 



- m 



[Ql] + _lim [xi] 

xi—>-0 



(D.43) 



The logarithmic singularity in the last line must again cancel for convergent integrals, and 
so we see that the only polynomials that appear in the limit are those in S(^x^y This finishes 
the proof of the second statement for the first integration. We stress that it is important 
that the integration region is [0, oo], because otherwise we have to take into account effects 
coming from the integration boundaries. 
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Let us now suppose that the two statements are true for the first r — 1 integrations. 
The set S(^xi,...,Xr-i) then consists of polynomials of the form = Q'^^Xr + R]^- Let us now 
compute the primitive with respect to Xr- We start by running the algorithm to bring 
the multiple polylogarithms in the integrand into canonical form. As this involves the 
application of the map (pxr ■> 

we obtain multiple polylogarithms of the form G(ai, ...,ayj;Xr) 

with 

ai&[^,-Rl/Ql} . (D.44) 

If we compute the primitive, we integrate over linear functions in from the set 5'(3,^ ^ , 
and it is easy to see that this does not change eq. (D.44). Using exactly the same argu- 
ment as for the first integration, we see that these multiple polylogarithms only contribute 
terms to the symbols whose entries are drawn from S(^xr,...xr-i)- addition we have mul- 
tiple polylogarithms that are independent of Xr, whose symbols involve those elements of 
S{x\,...xr-i) that are independent of Xr- As these functions do not change if we take the 
primitive with respect to x^, we see that they do not alter the conclusion. So the symbol 
of the primitive with respect to x^ must have all its entries drawn from S{^xi,...Xr-i)- Taking 
the limits Xr — )■ and x,- — )• oo just like for the first integration then finishes the proof. 

Having proved the two statements, we can show that our algorithm always terminates 
for the class of integral we consider. More precisely, we have to show that our algorithm 
can always produce the canonical form for the next integration step. This is done by 
applying the map cfy^^ , which requires all the entries in the symbol to be either independent 
of Xr or linear in Xr- By construction, this condition is always fulfilled for the elements of 
S{xi,...xr-i)- -'-t is easy to check that the same argument shows that the map (pxr+u which 
is called recursively by ^a;^ is well-defined, and so we can always find a canonical form for 
the integrand. 

D.4 A toy example 

In this section we discuss an example that illustrates how the algorithm described in the 
previous subsection can be used to compute parametric integrals. In particular, we used 
this algorithm to compute the parametric integrals (8.68) and (8.69). While these integrals 
are a straightforward application of the algorithm, intermediate expression are rather long, 
so we prefer to use a simpler integral where we can explicitly show all the steps. The 
integral we are going to consider is 

roo 

X(e) = / dxidx2dx3 x\ (1 + xif'-'^ ^2 ' (1 + xs)"^'"^ xf (1 + x^)-'-^ 

Jo (D.45) 

X (H- X2 -I- a;3 -I- xix^y'^'^'^ . 

The integral is finite as e ^ 0, and we want to compute the first few terms in the Taylor 
expansion 

X(e) = Jo + 2:1 e + X2 + 2:3 + O(e^) . (D.46) 
The first coefficient is trivial to compute, 

TT^ 2 

Z. = y-5. (D.47) 
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We will now illustrate our algorithm in detail on the coefficient Ii. If we integrate out the 
integration variables in the order (xi,X2,X3), we obtain 

S = {1 + Xi, 1 + X2, 1 + X3, 1 + X2 + X3 + X1X3} , 
S(x,) = {1 + X2, 1 + X3, 1 + X2 + X3} , (D.48) 

We see that all the sets are linear in all integration variables, so we perform the integrations 
one after the other. 

The coefficient Xi is given by the integral 

°° dxidx2dx3 



(1 + Xl)^ (1 + X2)^ (1 + X3) (1 + X2 + X3 + X1X3) 



3G (-1; xi) - 4G (-1; X2) - 3G (-1; X3) + G (0; xi) - G (0; X2) (d.49) 
+ 2G (0; X3) - 2G (-X3 - 1; X2) - 2G ("z^l^^inl; 



X3 

where we have already written all logarithms in terms of multiple polylogarithms, e.g., 

(X2 \ / X\ X3 
1 + + log 1 + — ■ 

1 + X3 y V 1 + X2 + X3 

= G(-l; X3) + G(-l - X3; X2) + G fz^^^^^; XI 

(D.50) 

It is easy to compute a primitive with respect to xi for the integrand of Xi, e.g.. 



/ 



-G(-l;x.)--Gf '"''"^'\ -l;x.). (D.51) 



1 + X2 + X3 + X1X3 X3 V X3 



In the following we only concentrate on this single term (which is in fact the most compli- 
cated one) to illustrate the procedure. All other terms can be dealt with in a similar way. 
Before we compute the limit of eq. (D.51) as xi 0, 00, let us comment about the symbol 
of the primitive. We have 



-X2 - X3 - 1 

G ,-l;xi 

X3 



(1 + Xi)0(l + X2) 



+ (1 + xi) (g) (1 + X2 + X3 + X1X3) - (1 + X2 + X3) (g) (1 + X2) (D.52) 

+ (1 + X2 + X3) (g X3 + (1+ X2 + X3 + X1X3) (8) (1 + X2) 
— (1 + X2 + X3 + X1X3) (g) X3 . 

We see that the entries in the symbol arc drawn from the set S U {xi, X2, X3}, as expected. 
The same is true for all other terms in the primitive. 
The limit xi is trivial, 

lim G ( -1; xi^ = . (D.53) 

xi-J-O V X3 / 
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The limit — > oo is obtained by letting xi = 1/xi and deriving the inversion relation for 
this multiple polylogarithm. This is equivalent to bringing G (—1, (— X2 — X3 — 1) / X3; 
into canonical form with respect to the ordering of variables (xi, X2, X3). We start by using 
our algorithm to construct a function G{xi,X2,X3) which is in canonical form and has the 
symbol given in eq. (D.52). We find 

g{xi,X2,X3) =G[ , ,-l;xi] -G( , , -, ,0;^i 

\ X2 + XS + I J \ X2 + X3 + I 

- G{0, -1; xi) + G(0, 0;xi)-G (-1, -X3 - 1; X2) ^^'^'^^ 

+ G (-1; X2) [G (0; X3) - G (-1; xg)] + G (0, -1; xg) - G (0, 0; xg) . 

It is easy to check that the symbol of G{xi,X2,xs) agrees with eq. (D.52). It would however 
be wrong to conclude that the functions are equal - they might differ by a rational multiple 
of ^2- Evaluating the difference numerically at a single point, we obtain 

G ( — — ^— ^ ) -e?(xi,X2,X3) = -1.64493406684822643... = -(2- (D.55) 



X3 xi^ 
Taking the limit xi — >■ 00 is now trivial, and we obtain 

G (-1, -^2-X3-l .^ \ ^ G(0,0;xi) - G(-l;x2)G(-l;x3) 



X3 

+ G (-1; X2) G (0; xg) - G (-1, -X3 - 1; X2) + G (0, -1; X3) - G (0, 0; X3) 
-C2 + 0(l/xi). 



(D.56) 



Note that the function has a logarithmic singularity for xi 00, which will cancel against 
similar contributions from other terms. Furthermore, note that the symbols of all the 
(finite) terms in eq. (D.56) have entries drawn form the set S'(j.j) U {x2,X3}, as expected. 
The previous steps can easily be implemented into a computer code and repeated for all 
the terms appearing in the primitive with respect to xi. 

The result of the integration over xi is, by construction, already in canonical form with 
respect to (x2,X3), and we can immediately compute the primitive with respect to X2, e.g., 



/ 



"^"^^ G(-l;x2) = G(-l-X3,-l;x2). (D.57) 



1 + X2 + X3 



The limit X2 — >■ is again trivial, while the limit X2 — )• 00 can again be computed by letting 
X2 = l/x2 and constructing a function ^(x2,X3) in canonical form with the same symbol. 
We find 



gi^.,X3) =G^-^,-l;.2j -G^--^,0;.2j -G(0,-l;.2) 
+ G(0,0;x2)-G(0,-l;x3) . 

Numerical evaluation at a single point immediately shows that we do not need to add any 
term proportional to (2- Taking the limit is now trivial, and we get 

G(-l - X3, -1; X2) = G(0, 0; X2) - G (0, -1; X3) + 0(l/x2) . (D.59) 
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We are finally only left with the integral over X3. The primitive involves integrals like 

dX3 



I 



I + X3 



G(-l,0;a;3) = G(-l,-l,0;x3). (D.60) 



The limit — > is trivial and for the limit X3 — )■ 00 we construct a function G{x2) with 
the same symbol as G(— 1, —1, 0; I/X3). We find 

g{x3) = G(-l, 0, 0; X3) - G{-1, -1, 0; X3) + G{0, -1, 0; X3) - G{0, 0, 0; X3) • (D.61) 

Next we have to determine terms proportional to (2- This is achieved by acting with A2,i 
on the difference, 



A2,i[G(-l,-l,0;l/x3)-g(x3)] 

= [G (-1, 0; X3) + G(-l, 0; X3) - G(0, 0; X3)] ® X3) 



+ 



G(0, 0; X3)-G{ -1, 0; ^ ) - G(-l, 0; ^3) 



(D.62) 



G(0;x3). 



Evaluating the first entries numerically at a single point reveals 



(D.63) 



A2,i -1, 0; I/X3) - g{x3)] = -C2 <8 X3) + C2 ® G(0; X3) 

= A2,i[-C2G(-l;x3) + C2G(0;x3)] . 

Finally, evaluating the full difference at a single point, we obtain 
G(-l, -1, 0; I/X3) - [g{x3) - C2 X3) + C2 G(0; X3)] = 1.20205690315959 ... 

(D.64) 

= Ca • 

Taking the limit X3 is now trivial, and we finally get 

277-2 5 

2:i = -5C3 + -^ + 3. (D.65) 

The higher terms in the e expansion can be obtained in exactly the same way. For this 
particular integral we find for example 

1497r^ 167r2 157 



216 9 6 /j3ggx 

^ __910 1497r^ 607 2777r^ 297r^ 1175 ^ ' ' 
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